Efficient convolution in an environment that enforces tiles

ABSTRACT

A method comprising: receiving an input tensor having a shape defined by [n1, . . . , nk], where k is equal to a number of dimensions that characterize the input tensor; receiving tile tensor metadata comprising: a tile tensor shape defined by [t1, . . . , tk], and information indicative of an interleaving stride to be applied with respect to each dimension of the tile tensor; constructing an output tensor comprising a plurality of the tile tensors, by applying a packing algorithm which maps each element of the input tensor to at least one slot location of one of the plurality of tile tensors, based on the tile tensor shape and the interleaving stride, wherein the interleaving stride results in non-contiguous mapping of the elements of the input tensor, such that each of the tile tensors includes a subset of the elements of the input tensor which are spaced within the input tensor according to the interleaving stride.

BACKGROUND

The invention relates generally to the field of encryption/decryptionschemes, algorithms, techniques, methods, computer programs, andsystems.

Fully homomorphic encryption (FHE) provides malleable ciphertext, whichallows operations to be performed on encrypted data, without firstdecrypting it. For example, if d denotes some data, E(d) denotes theencryption of d using a FHE scheme. Given E(d) and the scheme's publickey, it is possible to compute E(ƒ(d)) for any function ƒ, withoutknowing the decryption key, and without learning anything about d. Theresulting computations are left in an encrypted form which, whendecrypted, results in an identical output to that produced had theoperations been performed on the unencrypted data.

One potential use for FHE is in the case of outsourcing storage andprocessing of sensitive data, while preserving privacy. For example,using an outsourced service to compute a classification prediction overmedical data (e.g., medical images), while preserving the privacy of thepatient. Using an FHE scheme allows a data owner to send their data inan encrypted form to a cloud service that hosts a trained classifier.The encryption ensures that the data remains confidential, since thecloud service will not have access to the private key needed to decryptthe data. The cloud service will be capable of applying the trainedclassifier to the encrypted data to make encrypted predictions, andreturn the predictions in encrypted form to the data owner.

However, running large neural networks using FHE only is stillconsidered a computationally expensive task. This barrier forces usersto search for other secure alternatives instead of enjoying theadvantage of solutions that rely only on FHE.

The foregoing examples of the related art and limitations relatedtherewith are intended to be illustrative and not exclusive. Otherlimitations of the related art will become apparent to those of skill inthe art upon a reading of the specification and a study of the figures.

SUMMARY

The following embodiments and aspects thereof are described andillustrated in conjunction with systems, tools and methods which aremeant to be exemplary and illustrative, not limiting in scope.

There is provided, in an embodiment, a system comprising at least onehardware processor; and a non-transitory computer-readable storagemedium having stored thereon program instructions, the programinstructions executable by the at least one hardware processor to:receive an input tensor, wherein the input tensor has a shape defined by[n₁, . . . , n_(k)], where k is equal to a number of dimensions thatcharacterize the input tensor, receive tile tensor metadata comprisingat least: (i) a tile tensor shape defined by [t₁, . . . , t_(k)], and(ii) information indicative of an interleaving stride to be applied withrespect to each dimension of the tile tensor, and construct an outputtensor comprising a plurality of the tile tensors, by applying a packingalgorithm which maps each element of the input tensor to at least oneslot location of one of the plurality of tile tensors, based, at leastin part, on the tile tensor shape and the interleaving stride, whereinthe interleaving stride results in non-contiguous mapping of theelements of the input tensor, such that each of the tile tensorsincludes a subset of the elements of the input tensor which are spacedwithin the input tensor according to the interleaving stride.

There is also provided, in an embodiment, a computer-implemented methodcomprising: receiving an input tensor, wherein the input tensor has ashape defined by [n₁, . . . , n_(k)], where k is equal to a number ofdimensions that characterize the input tensor; receiving tile tensormetadata comprising at least: (i) a tile tensor shape defined by [t₁, .. . , t_(k)], and (ii) information indicative of an interleaving strideto be applied with respect to each dimension of the tile tensor;constructing an output tensor comprising a plurality of the tiletensors, by applying a packing algorithm which maps each element of theinput tensor to at least one slot location of one of the plurality oftile tensors, based, at least in part, on the tile tensor shape and theinterleaving stride, wherein the interleaving stride results innon-contiguous mapping of the elements of the input tensor, such thateach of the tile tensors includes a subset of the elements of the inputtensor which are spaced within the input tensor according to theinterleaving stride.

There is further provided, in an embodiment, a computer program productcomprising a non-transitory computer-readable storage medium havingprogram instructions embodied therewith, the program instructionsexecutable by at least one hardware processor to: receive an inputtensor, wherein the input tensor has a shape defined by [n₁, . . . ,n_(k)], where k is equal to a number of dimensions that characterize theinput tensor; receive tile tensor metadata comprising at least: (i) atile tensor shape defined by [t₁, . . . , t_(k)], and (ii) informationindicative of an interleaving stride to be applied with respect to eachdimension of the tile tensor; and construct an output tensor comprisinga plurality of the tile tensors, by applying a packing algorithm whichmaps each element of the input tensor to at least one slot location ofone of the plurality of tile tensors, based, at least in part, on thetile tensor shape and the interleaving stride, wherein the interleavingstride results in non-contiguous mapping of the elements of the inputtensor, such that each of the tile tensors includes a subset of theelements of the input tensor which are spaced within the input tensoraccording to the interleaving stride.

In some embodiments, the program instructions are further executable tostore, and the method further comprises storing the output tensor andthe tile tensor metadata.

In some embodiments, the program instructions are further executable tounpack, and the method further comprises unpacking, the input tensorfrom the stored output tensor, based on the tile tensor metadata.

In some embodiments, the tile tensor metadata further comprises areplication parameter, and wherein the packing algorithm is configuredto perform, based on the replication parameters, replication of theelements of the input tensor, such that each of the elements of theinput tensor is mapped to multiple slot locations along a dimension ofone of the tile tensors.

In some embodiments, the program instructions are further executable to,and the method further comprises: receive a filter associated with aconvolution computation over the input tensor; compute the convolutionby applying a multiplying operator which multiplies the filterelement-wise over each of the tile tensors in the output tensor; apply asummation algorithm that sums over the results of the multiplying; andoutput a result of the applying of the summation algorithm as a resultof the convolution.

In some embodiments, the convolution is part of a neural networkinference.

In some embodiments, the tensor tiles are homomorphic encryptionciphertexts.

In addition to the exemplary aspects and embodiments described above,further aspects and embodiments will become apparent by reference to thefigures and by study of the following detailed description.

BRIEF DESCRIPTION OF THE FIGURES

Exemplary embodiments are illustrated in referenced figures. Dimensionsof components and features shown in the figures are generally chosen forconvenience and clarity of presentation and are not necessarily shown toscale. The figures are listed below.

FIG. 1A is a block diagram of an exemplary system for data packingstructure for fully homomorphic encryption (FHE) schemes, according tosome embodiments of the present disclosure;

FIG. 1B presents a schematic of the packing optimizer, according to someembodiments of the present disclosure;

FIG. 1C is a flowchart of the functional steps in an example methodaccording to some embodiments of the present invention;

FIG. 2 is a diagram of a basic convolutional neural networkarchitecture.

FIGS. 3A-3C show tile tensor packing operations, according to someembodiments of the present disclosure;

FIG. 4A illustrates an interleaved tiling process, according to someembodiments of the present disclosure;

FIGS. 4B and 4C show two 8-parallel kernel evaluations using a 2×2filter or kernel performed over a matrix, according to some embodimentsof the present disclosure;

FIGS. 5A-5C illustrate an interleaved tiling operation, according tosome embodiments of the present disclosure;

FIGS. 6A-6D are another illustration of an interleaved tiling process ofthe present disclosure, according to some embodiments of the presentdisclosure; and

FIG. 7 is an illustration of experimental results with respect toCryptoNets implementation using tile tensors.

DETAILED DESCRIPTION

Disclosed herein is a technique, embodied in a method, system and acomputer program product, for data packing in the context of fullyhomomorphic encryption (FHE) schemes, which permits efficient high-leveltensor manipulation operations, such as convolution, over the encrypteddata, while reducing computational overhead.

By way of background, an FHE scheme is an encryption scheme that allowsto evaluate any circuit, and in particular any function, on encrypteddata. The FHE scheme receives as input data a vector M[s] and returns(ENC) a ciphertext. The FHE also generates a secret encryptionkey/public encryption key pair associated with the encrypted input data.The created cyphertext has a number of slots s that is determined duringthe key generation. The cyphertext may be deciphered (DEC) using thegenerated private encryption key, to return an s-dimensional vector,wherein M=Dec(Enc(M)). The functions for addition (Add), multiplication(Mul) and rotation (Rot) are then defined as:Dec(Add(Enc(M),Enc(M′)))=M+M′Dec(Mul(Enc(M),Enc(M′)))=M*M′Dec(Rot(Enc(M),n))(i)=M((i+n)mods)

Some FHE schemes, such as CKKS (see J. Cheon, et al., “Homomorphicencryption for arithmetic of approximate numbers,” in Proceedings ofAdvances in Cryptology—ASIACRYPT 2017. Springer Cham, 11 2017, pp.409-437), operate on ciphertexts in a homomorphic SIMD fashion. Thismeans that a single ciphertext encrypts a fixed size vector, and thehomomorphic operations on the ciphertext are performed slot-wise on theelements of the plaintext vector. To utilize the SIMD feature, more thanone input element must be packed and encrypted in each ciphertext. Thepacking method can dramatically affect the latency (i.e., time toperform computation), throughput (i.e., number of computations performedin a unit of time), communication costs, and memory requirements.However, deciding which packing to use is difficult and requires expertknowledge, and the more efficient packing may not be the trivial one.Moreover, different packing schemes offer different tradeoffs in termsof optimization goals. As the size of the FHE code increases, it becomesharder to find the optimal packing. For example, finding the bestpacking for a large neural network inference algorithm is a challengingtask, because the input is typically a four- or five-dimensional tensor,and the computation involves a long sequence of operations such asmatrix multiplication and convolution.

Accordingly, in some embodiments, the present disclosure provides for anFHE data packing technique which provides for high-level tensormanipulation operations, such as convolution.

In some embodiments, the present technique uses a data packing structuretermed ‘tile tensors.’ A tile tensor allows users to store tensors ofarbitrary shapes and sizes. The tile tensor automatically packs thetensor data into a collection of vectors of fixed size, as required inFHE environments, using a wide variety of configurable options. In someembodiments, tile tensors of the present disclosure also offer a set ofoperators to manipulate the tensor in its packed form. In someembodiments, tile tensors of the present disclosure supports the usageof operators on encrypted data, wherein the operators may be implementedusing generic algorithms that can work with any packing arrangement,agnostically of the packing arrangement selected internally. Thus, thepresent disclosure provides for a packing-obliviousprogramming-framework that allows users to concentrate on the design ofthe algorithms to be inferenced on the data, rather than data packingdecisions.

In some embodiments, the present disclosure further provides for apacking optimizer which operates in conjunction with the tile tensordata structure of the present disclosure. In some embodiments, thepresent optimizer searches for the optimum configuration for the tiletensor, given the user requirements and preferences. In someembodiments, the optimizer estimates the time and memory needed to run agiven function for every option, and returns the configuration whichoptimizes a given objective, whether latency, throughput, or memory. Insome embodiments, the present optimizer can be used to improve latencyfor small networks, adapt to various batch sizes, and scale up to muchlarger networks.

In some embodiments, the present disclosure provides a solution forconvolution computations, which are a popular building block in thefield of neural networks and machine learning models.

In some embodiments, the present disclosure is particularly useful inthe context of image processing. For example, in some embodiments, thepresent disclosure supports convolution over unlimited image size (e.g.,the total number of pixels) and reduces the number of rotations neededin the order of filterSize*sqrt (#imageSize−#slots). In someembodiments, the present disclosure can handle images with multiplechannels, convolution of kernels or filters, and batch processing ofmultiple images, and has easily tuned parameters to optimize forspecific scenarios.

In some embodiments, the present disclosure provides for a packingscheme which uses interleaved tiling, wherein input tensor values thatare adjacent to each other are stored within the same slot of differenttile. This way, when one wishes to calculate the convolution of a smallfilter over a large matrix, the calculations for the area covered by thefilter are much simpler, because all values are in different tiles butare within the same slot.

In some embodiments, the present disclosure interprets a ciphertext as amatrix of size T_(x) by T_(y), referred to as a ‘tile.’ If the inputtensor dimensions are I_(x) by I_(y), then a grid of

${Ex} = {{{{ceil}\left( \frac{Ix}{Tx} \right)}{}{by}{Ey}} = {{ceil}\left( \frac{Iy}{Ty} \right)}}$tiles are required to store the entire image. Input tensor values thatare adjacent or in close proximity within the same region may be storedsuch that they map to the same slot of different tiles. Specifically,the image pixel at (i,j) may be mapped to tile at position (i % Ex,j %Ey), and to the slot at position

$\left( {{{floor}\left( \frac{i}{Tx} \right)},{{floor}\left( \frac{j}{Ty} \right)}} \right).$

When calculating the convolution of a small filter over a large matrix,then all the calculations for the area covered by the filter are withinthe same slot. When calculations require crossing regions, i.e., mixingdifferent slot positions, the tiles may be rotated so that again all ofthe values are aligned within the same slot position. The exact numberof rotations can be calculated using the following method. Given animage size I_(x), I_(y), a filter having dimensions F_(x), F_(y), and atile sizes T_(x), T_(y), the number of rotations will be

${\left( {{Fx} - 1} \right)*\frac{Iy}{Ty}} + {\left( {{Fy} - 1} \right)*\frac{Ix}{Tx}} + {\left( {{Fx} - 1} \right)*{\left( {{Fy} - 1} \right).}}$Thus, the number of needed rotations is the square root of the size ofthe image, rather than of the same order as the size of the image. Thepresent method, combined with tile tensor data structures, can beextended with additional dimensions, thus it also allows handling ofmultiple images (batch processing), multiple image channels, andmultiple filters more efficiently. For example, sub-tile tensors may beused, where C is the number of channels, and B is the number of batches,wherein a sub-tile tensor contains the tensor [C,*, B], where *indicates a replicated dimension. The corresponding sub tile tensorshape would be

$\left\lbrack {\frac{C}{t1},{\frac{*}{t2}{,\frac{B}{t3}}}} \right\rbrack$for some tile dimensions [t1, t2, t3], and the full tile tensor,including all pixels, would be

$\left\lbrack {\frac{C}{t1},{\frac{*}{t2}{,\frac{B}{t3},\left. \frac{Ix}{t4} \right.\sim,\left. \frac{Iy}{t5} \right.\sim}}} \right\rbrack,$where the ˜ indicate interleaved dimensions. A filter's pixel is asub-tile tensor containing the tensor [C, F,*], with F denoting the sizeof the filter. Multiplying and summing over the C dimension results in atensor [*, F, B]. The full tile tensor shape for all filter pixels is

$\left\lbrack {{Fx},{Fy},\frac{C}{t1},\frac{F}{t2},{\frac{*}{t3}{,{\frac{*}{t4}{,\frac{*}{t5}}}}}} \right\rbrack.$Thus even if the image does fit inside the same ciphertext, but with allthe channels and batches it does not fit inside a single ciphertext, thepresent disclosure is able to avoid duplications and reduce rotations.Hardware and Software Environment

FIG. 1A is a block diagram of an exemplary system 100 for data packingstructure for fully homomorphic encryption (FHE) schemes, according tosome embodiments of the present disclosure.

System 100 may include one or more hardware processor(s) 102, arandom-access memory (RAM) 104, and one or more non-transitorycomputer-readable storage device(s) 106. As shown in FIG. 1A, system 100is an embodiment of a hardware and software environment for use withvarious embodiments of the present invention.

Storage device(s) 106 may have stored thereon program instructionsand/or components configured to operate hardware processor(s) 102. Theprogram instructions may include one or more software modules, such as adata packing module 108 and/or an optimizer module 110. The softwarecomponents may include an operating system having various softwarecomponents and/or drivers for controlling and managing general systemtasks (e.g., memory management, storage device control, powermanagement, etc.), and facilitating communication between varioushardware and software components. System 100 may operate by loadinginstructions of data packing module 108 and/or optimizer module 110 intoRAM 104 as they are being executed by processor(s) 102.

FIG. 1B is a block diagram of optimizer module 110. In some embodiments,optimizer module 110 is responsible for finding the most efficientpacking arrangement for a given computation, as well as the optimalconfiguration of the underlying FHE library, e.g., in connection withthe application of neural network operation to the FHE data.

In some embodiments, optimizer module 110 of the present disclosurereceives as input the model architecture of a neural network to beapplied. Optimizer module 110 automatically converts it to an FHEcomputation with optimal packing and optimal FHE library configuration.Users can further define packing constraints, such as the requiredsecurity level or maximal memory usage, and choose an optimizationtarget, whether to optimize for CPU time, latency or throughput, oroptimize for memory usage.

In some embodiments, optimizer module 110 selects from among differentpossible configurations of the FHE library, as well as different packingtechniques to support certain operators. In some embodiments, optimizermodule 110 also selects the tile shape, i.e., the values of t₁, t₂, . .. , in the tile tensor shapes. For example, consider an FHE schemeconfigured to have 16,384 slots in each ciphertext. Assuming a desiredconvolution operator uses five-dimensional tiles, the number of possibletuples t₁, . . . , t₅ such that

${{\prod}_{i}t_{i}} = {{16,384{is}\begin{pmatrix}{{\log_{2}(16,384)} + 5 - 1} \\{5 - 1}\end{pmatrix}} = {306{0.}}}$

In some embodiments, optimizer module 110 comprises three main units:the configuration generator 110 a, the cost evaluator 110 b, and thesimulator 110 c. in some embodiments, the optimization process may beginby a user providing a file (e.g., a JSON file) that contains the detailsof the model architecture to be applied in runtime. The configurationgenerator 110 a generates a list of all possible configurations,including the packing details and FHE configuration details applicablefor this architecture. The simulator unit 110 b tests every suchconfiguration and outputs one or more of the following data for each:the computation time of the different stages including encrypting themodel and input samples, running inference, and decrypting the results;throughput; memory usage of the encrypted model; input; and output.Optimizer module 110 passes this data to the cost evaluator 110 c forevaluation. Finally, optimizer module 110 returns the configurationoption that yields the optimal cost to the user, together with thesimulation output profile.

In some embodiments, the configuration generator unit 110 a of optimizermodule 110 receives the model architecture to be applied in runtime, andgenerates all applicable configurations for it. The generator unit 110 awill then create multiple complete configurations by exploring allpossible tile shapes. The generator unit 110 a explores possible tileshape using one of two strategies. The first involves brute forcing overall valid options for tile shapes. Since these may be numerous, a secondstrategy searches using a “steepest ascent hill climbing” local searchalgorithm. The local search starts with a balanced tile shape, where thenumber of slots in every dimension is of the same order. This is aheuristic designed to avoid evaluating tile shapes that are likely to becomputationally costly at the beginning of the search. Thus, all theneighbor tile shapes of the current shape can be iteratively evaluated,wherein the best-improving neighbor is selected, as long as one exists.In some embodiments, two tile shapes may be considered as neighbors,where one shape may be obtained from the other by multiplying ordividing the size of some of its dimensions by two. A tile shape may bedeemed to be better over another tile shape based on the costs receivedfrom the cost evaluator. Using the local search algorithm highly speedsup the search process, and often results in a global optimum.

In some embodiments, the simulator unit 110 b receives as inputs themodel architecture to be applied in runtime, and a configuration optionfrom configuration generator 110 a. At this stage, the configuration maybe evaluated by running it on encrypted input under FHE. To reducecomputational costs, the simulator 110 b may use pre-calculatedbenchmark values such as the CPU time of every HE operation and thememory consumption of a tile (i.e., the memory consumption of a singleciphertext). Then, simulator 110 b evaluates the model to be applied onmockup tile tensor objects using these benchmarks. These mockup tiletensors contain only meta data and gather performance statistics. Usingthis approach, the simulator 110 b can simulate an inference operationseveral order-of-magnitudes faster than when running the complete modelon encrypted data.

In some embodiments, the cost evaluation unit 110 c evaluates thesimulator 110 b output data considering the constraints and optimizationtargets, which may be user-provided. After testing all possibleconfigurations, the highest scoring configuration(s) is sent back asoutput.

System 100, as described herein, is only an exemplary embodiment of thepresent invention, and in practice may be implemented in hardware only,software only, or a combination of both hardware and software. System100 may have more or fewer components and modules than shown, maycombine two or more of the components, or may have a differentconfiguration or arrangement of the components. System 100 may includeany additional component enabling it to function as an operable computersystem, such as a motherboard, data busses, power supply, a networkinterface card, a display, an input device (e.g., keyboard, pointingdevice, touch-sensitive display), etc. (not shown). Moreover, componentsof system 100 may be co-located or distributed, or the system may beconfigured to run as one or more cloud computing “instances,”“containers,” “virtual machines,” or other types of encapsulatedsoftware applications, as known in the art.

The programs described herein are identified based upon the applicationfor which they are implemented in a specific embodiment of theinvention. However, it should be appreciated that any particular programnomenclature herein is used merely for convenience, and thus theinvention should not be limited to use solely in any specificapplication identified and/or implied by such nomenclature.

The descriptions of the various embodiments of the present inventionhave been presented for purposes of illustration, but are not intendedto be exhaustive or limited to the embodiments disclosed. Manymodifications and variations will be apparent to those of ordinary skillin the art without departing from the scope and spirit of the describedembodiments. The terminology used herein was chosen to best explain theprinciples of the embodiments, the practical application or technicalimprovement over technologies found in the marketplace, or to enableothers of ordinary skill in the art to understand the embodimentsdisclosed herein.

As shown in FIG. 1A, system 100 is an environment in which an examplemethod according to some embodiments of the present invention can beperformed.

FIG. 1C is a flowchart of the functional steps in an example methodaccording to some embodiments of the present invention. The method ofthe flowchart is a method of an interleaved packing of a tile tensorfrom an input tensor. This method and associated software for packingthe tile tensor will now be discussed, over the course of the followingparagraphs, with reference to of FIGS. 1A-1C.

Processing begins at step 120, where system 100 receives the data of aninput tensor A1 data set with k dimensions of sizes n₁, n₂ . . . n_(k).A tensor is a multi-dimension array, with k dimensions and shape [n₁,n₂, . . . n_(k)]. For example, if the input tensor is a matrix with 4rows and 5 columns, then k=2, n₁=4, n₂=5. If the input tensor is avector of length 7, then k=1, n₁=7. If the input tensor is a 3D array ofdimensions 10 by 20 by 30, then k=3, n₁=10, n₂=20, n₃=30.

In step 122, system 100 further receives tile tensor metadata, includinga tile tensor shape and an interleaving stride. The tile tensor shapeconsists of the shape of the dimensions of the input tensor A1, namelyk, n₁, n₂, . . . n_(k), and the packing details t₁, t₂, . . . , t_(k).The tile tensor data structure includes copy of the tile tensor shape aswell.

The required number of tile tensors can be computed as follows. Let

${e_{i} = {{ceil}\left( \frac{n_{i}}{t_{i}} \right)}},$where cell means rounding up to the nearest integer from above. Therequired number of tiles is e₁*e₂* . . . *e_(k).

In step 124, the instructions of system 100 cause packing module 108 toapply a packing algorithm, which maps each element of the input tensorto at least one slot location in the tile tensors, in a non-contiguousmanner, based on the tile tensor shape and the interleaving stride. Insome embodiments, the interleaving stride results in non-contiguousmapping of the elements of said input tensor, such that each of the tiletensors includes a subset of the elements of the input tensor which arespaced within the input tensor according to the interleaving stride.

In step 126, the instructions of system 100 cause packing module 108 toconstruct an output tensor comprising a plurality of the tile tensors.

In step 128, the instructions of system, 100 may cause it to store,e.g., in storage device 106, the output tensor and the tile tensormetadata.

In step 130, the instructions of system, 100 may cause it to unpack theinput tensor from the stored output tensor, based on the tile tensormetadata.

Tensor Basic Operations

As used herein, the term ‘tensor’ is synonymous with multi-dimensionalarray, as this is common in the AI domain.

A k-dimensional tensor may be denoted by [n₁, n₂, . . . , n_(k)], where0≤n_(i) is the size of the i'th dimension. For example, the shape of a5×6 matrix M is [5,6]. For a tensor R[n₁, . . . , n_(k)], R(j₁, j₂, . .. , j_(k)) may be used to refer to a specific element, where0≤j_(i)<n_(i).

Matrix multiplication operation shall be denoted herein without amultiplication symbol, e.g., M₁M₂ stands for the product of M₁ and M₂.The transpose operation of a matrix M may be denoted by M^(T), and tags(e.g., M′,M″) may be used to denote different objects. The operationsM₁+M₂ and M₁*M₂ refer to element-wise addition and multiplication,respectively.

The tensors A[n₁, . . . , n_(k)] and B[m₁, . . . , m_(k)] are said tohave compatible shapes if m_(i)=n_(i) or either n_(i)=1 or m_(i)=1, fori≤k. Their mutual expanded shape is [max{n_(i), m_(i)}]_(i≤k).

When a tensor A has more dimensions than a tensor B, their dimensionscan be matched by expanding B with dimensions of size 1. This results inequivalent tensors up to transposition. For example, both tensors V[b]and V[b, 1] represent column vectors, while V[1, b]=V^(T) represents arow vector.

The broadcasting operation takes two tensors with compatible butdifferent shapes and expands every one of them to their mutual expandedshape. For a tensor A[n₁, . . . , n_(k)] and a tensor shape s=[m₁, . . ., m_(k)] with n_(i)∈{1, m_(i)} for each i=1, . . . , k, the operationC=broadcast(A, s) replicates the content of A along the r dimensionm_(r) times for every r=1, . . . , k and n_(r)=1<m_(r). The outputtensor C is of shape s.

The tensors A[3,4,1] and B[1,4,5] have compatible shapes. Their mutualexpanded shape is s=[3,4,5] and broadcast(A, s) has the same shape s asbroadcast(B, s).

Element-wise operations, such as addition (A+B) and multiplication (A*B)on two tensors with compatible shapes A, B, may be performed by firstusing broadcasting to expand them to their mutual expanded shape andthen performing the relevant element-wise operation. Thus, for a tensorA[n₁, . . . , n_(k)], the operation B=sum(A, t) sums the elements of Aalong the t-th dimension and the resulting tensor B has shape [n₁, . . ., n_(t-1), 1, . . . , n_(k)] and

${{B\left( {j_{1},\ldots,j_{t - 1},0,\ldots,j_{k}} \right)} = {{\sum}_{l = 0}^{n_{t} - 1}{A\left( {j_{1},\ldots,j_{t - 1},l,\ldots,j_{k}} \right)}}},$for all j_(i)<n_(i) for i∈{1, 2, . . . , k}\{t},

Using broadcasting and summation, common algebraic operators can beperformed. For two matrices M₁ [a, b], M₂ [b, c] and the column vectorV[b, 1], matrix-vector multiplication may be performed usingM₁V=sum(M₁*V^(T), 2), where M₁ and V^(T) have compatible shapes with themutual expanded shape of [a, b]. Matrix-matrix multiplication may beperformed using M₁M₂=sum(M₁, *M_(2′),2), where M_(1′)=M₁[a, b, 1] andM_(2′)=M₂[1, b, c].

Tile Tensors

In some embodiments, a ‘tile tensor’ is a data structure containing anexternal tensor as data and a tile tensor shape as meta data. An‘external tensor’ is a tensor in which each element is a tile.

In some embodiments, a tile tensor of the present disclosure is a datastructure that packs tensors in fixed size chunks, as required for FHE,and allows them to be manipulated similarly to regular, i.e.,non-encrypted, tensors. In other words, performing operations on thetile tensors of the present disclosure may be considered to beequivalent (or “homomorphic,” as is the mathematical term) to directlyoperating on the input tensors.

As the term is used herein, a “tensor” is any multi-dimensional array ofnumbers. One special case for tensors is when the array is onedimensional, in which case the array is generally referred to as avector. In the case of a two-dimensional vector, it is commonly referredto as a matrix. It is understood in the art that the numbers making up atensor can be integers or floating-point numbers or complex numbers.However, in some tensors, the elements of a vector may not take the formof numbers, but, rather, these elements may take other forms, such ascharacters, strings, or other objects. In the specific examplesdiscussed in this document, it will generally be assumed that thetensors contain elements in the form of numbers unless otherwise noted.

It is known that tensors can be packed into tiles and, further, thatmathematical operators can be performed on these tensors while in packedform. As used herein, a “tile” is any contiguous block of data storagein a data storage device (for example, a volatile memory type datastorage device) of fixed size, capable of holding n numbers, where n issome fixed size determined by system configuration. In this way, alltiles will have the same size of n numbers.

It is understood by those of skill in the art that mathematicaloperations can be performed on the elements of tiles (for example, thenumbers making up a tile). Such a mathematical operation is typicallyperformed element-wise on the respective numbers of the respectivetile(s) being subject to the mathematical operation. For example, iftile T₁ has the numbers (x₁, x₂, . . . x_(n)) and tile T₂ has thenumbers (y₁, y₂, . . . y_(n)), then T₁+T₂ is a tile containing (x₁+y₁,x₂+y₂, . . . x_(n)+y_(n)). And similarly, T₁*T₂ is a tile containing(x₁*y₁, x₂*y₂, . . . , x_(n)*y_(n)). Tiles can also be rotated by an anyoffset r, which means moving each element r slots to the left (or to theright, if r is negative), and first elements rotate back to being last.For example, if T₁=(1,2,3,4,5,6,7,8) then T_1 rotated by 2 is(3,4,5,6,7,8,1,2) and rotate by −2 is (7,8,1,2,3,4,5,6).

The process of creating a tile tensor from an input tensor is termedpacking. The inputs to this process are the tensor itself and a tiletensor shape, e.g., meta-data indicating the packing scheme. The outputis a tile tensor data structure containing a set of tiles filled withdata copied from the tensor and arranged within the tiles according tothe tile tensor shape. The tile tensor data structure includes a copy ofthe tile tensor shape as well. The tensor can be retrieved back from thetile tensor using a process called “unpacking.” The tile tensor shape isused to identify how the tensor's elements are arranged within thetiles, and copied back into a tensor. Given one or more tile tensors,operators can be applied to them. These operators may change both thecontent of the tiles, and the tile tensor shapes of the tile tensors.

In some embodiments, an ‘external tensor’ as used herein is ak-dimensional tensor wherein each of its elements is itself ak-dimensional tensor, all having an identical shape. These internaltensors are referred to as ‘tiles’ or ‘tile tensors,’ their shape is the‘tile shape,’ and the shape of the external tensor is the ‘externalshape.’ A slot in E is identified by E(a₁, . . . , a_(k))(b₁, . . . ,b_(k)) where a_(i) are the external indices of a tile, and b_(i) are theinternal indices inside the tile.

A k-dimensional ‘tile tensor shape’ is comprised of an external shape[e₁, . . . , e_(k)], tile shape [t₁, . . . t_(k)], original shape [n₁, .. . , n_(k)], replication counts [r₁, . . . r_(k)], interleaved Booleanindicator [l₁, . . . , l_(k)], and unknown Boolean indicators [u₁, . . ., u_(k)]. It is required that ∀_(i)(r_(i)=1∨n_(i)=1)∧(max(r_(i),n_(i))≤e_(i)t_(i)).

Given a tile tensor shape S, an external tensor E, a specific slot in Especified by external indices (a₁, a_(k)), and internal indices (b₁, . .. b_(k)), then this slot is associated with the logical indices (c₁, . .. , c_(k)) with respect to S, computed as follows: For i=1, . . . , k,if the interleaved indicator l_(i) is true, then c_(i)=b_(i)e_(i)+a_(i)else c_(i)=a_(i)t_(i)+b_(i).

A tile tensor shape S is valid for an external tensor E if theirexternal shapes and tile shapes match, and there exists a tensor T[n₁, .. . , n_(k)] such that for T₁=broadcast(T, [n₁n₂, . . . , n₂r₂]) itholds that E(a₁, . . . , a_(k))(b₁, . . . , b_(k))=T₁(c₁, . . . , c_(k))for all slots with internal, external, and logical indices a_(i), b_(i),c_(i), such that ∀_(i)c_(i)≤n_(i)r_(i). For all other slots of E, if∀_(i)((c_(i)≥r_(i)n_(i))→¬u_(i))) then these slots are set to zero. T isthe packed tensor.

Tile tensor is a pair (E, S) where E is an external tensor and S a tiletensor shape that is valid for it.

Given a tile tensor T_(A)=(E, S) the operator unpack(E) results with thepacked tensor of T_(A).

Given a tensor A and a tile tensor shape S whose original shape matchesthe shape of A, then the pack operator pack(A, S) results with a tiletensor T_(A)=(E, S) such that A is the packed tensor of T_(A).

A tile tensor shape can be specified with a special notation involving alist of symbols. Each element in the list specifies the details of onedimension.

$\frac{n_{i}}{t_{i}}$specifies the original and tile shape along this dimension, and r_(i)=1,

${e_{i} = \frac{n_{i}}{t_{i}}},$l_(i)=u_(i)=false.

$\frac{*r_{i}}{t_{i}}$further specifies the replication count and n_(i)=1, and

$\frac{*}{t_{i}}$specifies n_(i)=1, r_(i)=t_(i).

$\frac{n_{i} \sim}{t_{i}}$specifies l_(i)=true, and

$\frac{n_{i} \sim e_{i}}{t_{i}}$specifies a value for e_(i) other than the default mentioned above. Forany of the above mentioned options a “?” symbol above the line indicatesu_(i)=true.Tile Tensor Data Structure

In some embodiments, a basic tiling process according to someembodiments of the present disclosure may comprise taking a tensor A[n₁,n₂, . . . , n_(k)], and breaking it down into equal-size blocks, ortiles, each having the shape [t₁, t₂, . . . , t_(k)].

FIG. 3A shows a tile tensor packing operation. Input tensor 302, e.g., aone-dimensional vector comprising 8 values 00-07, may be reinterpretedas a multi-dimensional array 304, e.g., [2,4].

In some embodiments, a tensor E[e₁, e₂, . . . , e_(k)] may then beconstructed, which may be termed an ‘external tensor,’ such that eachelement of E is a tile, and

${e_{i} = \frac{n_{i}}{t_{i}}}.$Thus, T=E(a₁, a₂, . . . , a_(k)) for 0≤a_(i)<e_(i), is a specific tilein E, and T(b₁, b₂, . . . , b_(k)) for 0≤b_(i)<t_(i) is a specific slotinside this tile. An element of the original tensor A(c₁, c₂, . . . ,c_(k)) will be mapped to tile indices

${a_{i} = \frac{c_{i}}{t_{i}}},$and indices inside the tile b_(i)=c_(i) mod t_(i). All other slots in Ethat were not mapped to any element of A will be set to 0. FIG. 3Bdemonstrates three examples for tiling a matrix M[5,6] using 8-slot tiletensors of shape [2,4]. As can be seen, the resulting shape of theexternal tensor is [3,2], and the tile shape is [2,4]. FIG. 3C showstiling a matrix M[9,10] using tiles of shape [2,4].

Tile tensor shapes may be denoted using the following notation. Forexample,

$\left\lbrack {\frac{n_{1}}{t_{1}},\frac{n_{2}}{t_{2}},\ldots,\frac{n_{k}}{t_{k}}} \right\rbrack$is a tile tensor shape specifying an input tensor of shape [n₁, . . . ,n_(k)], packed and tiled using tiles of shape [t₁, . . . , t_(k)]. Inthis notation, if t_(i)=1, then it can be omitted. For example,

$\left\lbrack {\frac{5}{1},\frac{6}{8}} \right\rbrack$can be written

$\left\lbrack {5,\frac{6}{8}} \right\rbrack.$

A tile tensor can be created using a pack operation that receives atensor A[n₁, . . . , n_(k)] to be packed and the desired tile tensorshape:

$T_{A} = {{{pack}\left( {A,\ \left\lbrack {\frac{n_{1}}{t_{1}},\ldots,\frac{n_{k}}{t_{k}}} \right\rbrack} \right)}.}$The pack operator computes the external tensor using the tiling processdescribed above, and stores along-side it the tile tensor shape, to formthe full tile tensor T_(A). A may be retrieved back using the unpackoperation: A=unpack (T_(A)). As with regular tensors, a tile tensorT_(A) may be denoted together with its shape:

${T_{A}\left\lbrack {\frac{n_{1}}{t_{1}},\ldots,\frac{n_{k}}{t_{k}}} \right\rbrack}.$

In some embodiments, tile tensors can be manipulated by operators bychanging the content of the multi-dimensional array of tiles and theaccompanying packing details and other meta data. The change is made insuch a way that the operator is equivalent to applying such an operatordirectly tensors they contain packed inside. For example, if tensor A₁is packed inside tile tensor TA₁, and tensor A₂ is packed inside tiletensor TA₁, then software instructions can apply the “add” operator onTA₁ and TA₂, obtaining a new tile tensor TA₃. Unpacking TA₃ will resultin a tensor A3, which is equal to A₁+A₂.

Replication

For some computations it is useful to have the tensor data replicatedseveral times inside the tile slots. The tile tensor shape indicatesthis by using the

$\frac{*}{t_{i}}$notation. It implies that n_(i)=1, but each element of the originaltensor is replicated t_(i) times along the i'th dimension. When packinga tensor A[n₁, . . . , n_(k)] and n_(i)=1, and with a tile tensor shapespecifying

$\frac{*}{t_{i}},$then the packing operation performs broadcast(A, [n₁, . . . , t_(i), . .. , n_(k)]) and tiles the result. The unpacking process shrinks thetensor back to its original size. The replications can either beignored, or an average of them can be taken; this is useful if the datais stored in a noisy storage medium, as in approximate FHE schemes.Unknown Values

When tensors are packed into tile tensors, unused slots are filled withzeroes, as shown, e.g., in FIG. 3B. When operators are applied on tiletensors, unused slots might get filled with arbitrary values. Althoughthese unused slots are ignored when the tile tensor is unpacked, thepresence of arbitrary values in them can still impact the validity orperformance of applying additional operators. To reflect this state, thetile tensor shape contains an additional flag per dimension, denoted bythe symbol “?”, indicating the presence of unknown values.

FIG. 3D shows packing tensor V[5, 1] into tile tensors using differenttile tensor shapes:

-   -   Panel A uses a tile tensor with the shape

$T_{V\lbrack{\frac{5}{2},\frac{1}{4}}\rbrack}.$

-   -   Panel B demonstrates tiling with replication, where the packing        process computes V′=broadcast(V, [5,4]) and uses tiles of shape

$T_{V\lbrack{\frac{5}{2},\frac{*}{4}}\rbrack}^{\prime}.$

-   -   Panel C uses a tile tensor with the shape

$T_{V\lbrack{\frac{5}{2},\frac{1?}{4}}\rbrack}^{''}.$

-   -    The “?” in the second dimension indicates that whenever the        valid range of the packed tensor is exceeded along this        dimension, there may be encountered arbitrary unknown values.        However, it still holds that V=unpack(T_(V)), as these unused        slots are ignored.        Operators

Tile tensor operators are homomorphic operations between tile tensorsand the packed tensors they contain. For two tile tensors T_(A) andT_(B), and a binary operator ⊙, it holds thatunpack(T_(A)⊙T_(B))=unpack(T_(A))⊙unpack(T_(B)). Unary operators aresimilarly defined.

Binary elementwise operators are implemented by applying the operator onthe external tensors tile-wise, and the tile tensor shape is updated toreflect the shape of the result. If the inputs have identical shape,then so do the results, e.g.,

$T_{M^{''}}\left\lbrack {\frac{5}{2},\frac{6}{4}} \right\rbrack$in FIG. 3C can be multiplied with an identically packed matrix

${T_{N}\left\lbrack {\frac{5}{2},\frac{6}{4}} \right\rbrack},$resulting in

${T_{R}\left\lbrack {\frac{5}{2},\frac{6}{4}} \right\rbrack},$where R=M*N. As with regular tensors, the tile tensor shapes need not beidentical, but compatible. Compatible tile tensor shapes have the samenumber of dimensions, and for each dimension specification they areeither identical, or one is

$\frac{*}{t_{i}}$and the other is

$\frac{n_{i}}{t_{i}}.$The intuition is that if the tensor is already broadcast inside thetile, it can be further broadcast to match any size by replicating thetile itself. For example, for

$T_{V^{\prime}}\left\lbrack {\frac{5}{2},\frac{*}{4}} \right\rbrack$panel B, TM_(″)*T_(V′) may be computed, resulting in

$T_{R^{\prime}} = {\left\lbrack {\frac{5}{2},\frac{6}{4}} \right\rbrack.}$T_(M″)+T_(V′), may also be computed, but this results in

${T_{R^{''}}\left\lbrack {\frac{5}{2},\frac{6?}{4}} \right\rbrack},$i.e., with unknown values in unused slots along the second dimension.This occurs because in T_(V), this dimension is filled with replicatedvalues, and after the addition they fill the unused slots of the result.Computing T_(M″)*T_(V) is illegal because their shapes are notcompatible.

The sum operator is also defined homomorphically: unpack(sum(T_(A),i))=sum(unpack(T_(A)), i). It works by summing over the external tensoralong the i'th dimension, then by summing inside each tile along thei'th dimension. In an FHE environment, the latter summation requiresusing the rotate-and-sum algorithm. Generally, the sum operator reducesthe i'th dimension and the resulting tile tensor shape changes to

$\frac{1?}{t_{i}}.$However, mere are some useful special cases. If t_(i)=1, then it isreduced to 1/1 or simply 1. When i is the smallest i such that t_(i)>1,the dimension reduces to

$\frac{*}{t_{i}},$i.e., the sum results are replicated. This is due to properties of therotate-and-sum algorithms. It is a useful property, since thisreplication is sometimes needed for compatibility with another tiletensor. For example, let T_(A) be a tile tensor with the shape

$\left\lbrack {4,\frac{3}{8},\frac{5}{16}} \right\rbrack.$Then sum(T_(A), 1) is of shape

$\left\lbrack {1,\frac{3}{8},\frac{5}{16}} \right\rbrack;$sum(T_(A), 2) is of shape

$\left\lbrack {4,{\frac{*}{8}{,\frac{5}{16}}}} \right\rbrack$and sum(T_(A), 3) is of shape

$\left\lbrack {4,\frac{3}{8},\frac{1?}{16}} \right\rbrack.$

Three other operators do not change the packed tensor, just the externaltensor and tile tensor shape. The clear(T_(A)) operator clears unknownvalues by multiplying with a mask containing ones for all used slots,i.e., it removes the “?” from the tile tensor shape. For example,

${{clear}\left( {T_{V^{''}}\left\lbrack {\frac{5}{2},\frac{1?}{4}} \right\rbrack} \right)} = {{T_{V}\left\lbrack {\frac{5}{2},\frac{1?}{4}} \right\rbrack}.}$The rep(T_(A), i) operator assumes the i'th dimension is

$\frac{1}{t_{i}},$and replicates it to

$\frac{*}{t_{i}},$using a rotate-and-sum algorithm. The flatten(T_(A), i, j) operatorflattens dimensions i through j assuming they are all replicated. Thisis done trivially by just changing the meta data, e.g.,

${flatten}\left( {{T_{A}\left\lbrack {\frac{3}{4},\frac{*}{8},\frac{*}{2},\frac{5}{32}} \right\rbrack},2,3} \right)$results with

$T_{A^{\prime}}\left\lbrack {\frac{3}{4},{\frac{*}{16}{,\frac{5}{32}}}} \right\rbrack$Higher Level Operators

Using elementwise operators and summation, various algebraic operationsmay be performed on tile tensors.

In some embodiments, given a matrix M[a, b] and a vector V[b], V may bereshaped into V[1, b] for compatibility, and pack both tensors into tiletensors as

${T_{M}\left\lbrack {\frac{a}{t_{1}},\frac{b}{t_{2}}} \right\rbrack},$and

${T_{V}\left\lbrack {\frac{*}{t_{1}},\frac{b}{t_{2}}} \right\rbrack},$for some chosen tile shape [t₁, t₂]. These may be then multiplied using:

$\begin{matrix}{{T_{R}\left\lbrack {\frac{a}{t_{1}},\frac{1?}{t_{2}}} \right\rbrack} = {{{sum}\left( {{{T\left\lbrack {\frac{a}{t_{1}},\frac{b}{t_{2}}} \right\rbrack}*{T_{V}\left\lbrack {\frac{*}{t_{1}}{,\frac{b}{t_{2}}}} \right\rbrack}},2} \right)}.}} & (1)\end{matrix}$

Eq(1) works for any value of a, b, t₁, t₂. This is because the tiletensor shapes of T_(M) and T_(V) are compatible, and therefore, due tothe homomorphism, this computes R[a, 1]=sum(M[a, b]*V[1, b], 2), whichproduces the correct result as explained above.

A second option is to initially transpose both M and V and pack them intile tensors

$T_{M}\left\lbrack {\frac{b}{t_{1}},\frac{a}{t_{2}}} \right\rbrack$and

${T_{V}\left\lbrack {\frac{b}{t_{1}},\frac{*}{t_{2}}} \right\rbrack}.$Now they may be multiplied as:

$\begin{matrix}{{T_{R}\left\lbrack {\frac{*}{t_{1}}{,\frac{a}{t_{2}}}} \right\rbrack} = {{{sum}\left( {{{T_{M}\left\lbrack {\frac{b}{t_{1}},\frac{a}{t_{2}}} \right\rbrack}*{T_{V}\left\lbrack {\frac{b}{t_{1}},\ \frac{*}{t_{2}}} \right\rbrack}},1} \right)}.}} & (2)\end{matrix}$

This computes the correct result using the same reasoning as before. Thebenefit here is that the result

$T_{R}\left\lbrack {\frac{*}{t_{1}},\frac{a}{t_{2}}} \right\rbrack$is replicated along the first dimension due to the properties of the sumoperator. Thus, it is ready to play the role of T_(V) in Eq(1) above,and two matrix-vector multiplications consecutively may be performedwithout any processing in between. The output of Eq(1) above can beprocessed to fit as input for Eq(5) below, using rep(clean(T_(R)),2).

The above reasoning easily extends to matrix-matrix multiplication asfollows. Given matrices M₁[a, b] and M₂[b, c], their product may becomputed using either of the next two equation, where in the second oneM₁ is transposed prior to packing. As before, the result of the secondfits as input to the first.

$\begin{matrix}{{T_{R}\left\lbrack {\frac{a}{t_{1}},\frac{1?}{t_{2}},\frac{c}{t_{3}}} \right\rbrack} = {{{sum}\left( {{{T_{M_{1}}\left\lbrack {\frac{a}{t_{1}},\ \frac{b}{t_{2}}\ ,\frac{*}{t_{3}}} \right\rbrack}*{T_{M_{2}}\left\lbrack \ {\frac{*}{t_{1}}{,\frac{b}{t_{2}}\ ,\ \frac{c}{t_{3}}}} \right\rbrack}},2} \right)}.}} & (3)\end{matrix}$ $\begin{matrix}{{T_{R}\left\lbrack {\frac{*}{t_{1}}{,\frac{a}{t_{2}},\frac{c}{t_{3}}}} \right\rbrack} = {{{sum}\left( {{{T_{M_{1}}\left\lbrack {\frac{b}{t_{1}},\frac{a}{t_{2}},\ \frac{*}{t_{3}}} \right\rbrack}*{T_{M_{2}}\left\lbrack {\frac{b}{t_{1}}\ ,{\frac{*}{t_{2}}{,\frac{c}{t_{3}}}}} \right\rbrack}},1} \right)}.}} & (4)\end{matrix}$

The product R[100,60]=Π_(i=1) ⁴ M_(i) of the four matrices M₁[100,90],M₂[90,80], M₃[80,70], and M₄[70,60] is computed by packing the matricesin tile tensors

${T_{M_{1}}\left\lbrack {\frac{90}{t_{1}},\frac{100}{t_{2}},\frac{*}{t_{3}}} \right\rbrack},{T_{M_{2}}\left\lbrack {\frac{90}{t_{1}},\frac{80}{t_{2}},\frac{*}{t_{3}}} \right\rbrack},{T_{M_{3}}\left\lbrack {\frac{70}{t_{1}},\frac{80}{t_{2}},\frac{*}{t_{3}}} \right\rbrack},$${and}{T_{M_{4}}\left\lbrack {\frac{70}{t_{1}},{\frac{*}{t_{3}}\ {,\frac{60}{t_{3}}}}} \right\rbrack}$and computing

$\begin{matrix}{{T_{X_{1}}\left\lbrack {\frac{90}{t_{1}},\ \frac{1?}{t_{2}},\ \frac{60}{t_{3}}} \right\rbrack} = {{sum}\left( {{T_{M_{2}}*\left( {{sum}\left( {{T_{M_{3}}*T_{M_{4}}},1} \right)} \right)},2} \right)}} \\{{T_{R}\left\lbrack {\frac{*}{t_{1}}{,\frac{100}{t_{2}},\frac{60}{t_{3}}}} \right\rbrack} = {{sum}\left( {{T_{M_{1}}*\left( {{rep}\left( {{{clean}\left( T_{X_{1}} \right)},2} \right)} \right)},1} \right)}}\end{matrix}$Interleaved Tiling

In some embodiments, the present disclosure provides for a tilingprocess wherein the tiles do not cover area of the tensor in asequential, contiguous manner, but rather are spread over the tensorusing equal strides.

Another option for tiling is denoted by the symbol “˜” in the tiletensor shape. This symbol indicates that the tiles do not cover acontiguous block of the tensor, but are spread out in equal strides. Ifthe dimensions are interleaved, an element of the original tensor A(c₁,c₂, . . . , c_(k)) will be mapped to tile indices a_(i)=c_(i)mode_(i),and indices inside the tile

$b_{i} = \frac{c_{i}}{e_{i}}$(where e_(i) is the size of the external tensor).

FIG. 4A illustrates an interleaved tiling process of the presentdisclosure. As can be seen, input tensor 402 is a matrix M[6,6]. In someembodiments, input matrix 402 may be, e.g., an image, wherein eachelement in the matrix may be, e.g., an image pixel.

Input tensor 402 is packed using tile tensors 404 a-404 f of shape[2,4], notated

${T_{M}\left\lbrack {\frac{6 \sim}{2},\frac{6 \sim}{4}} \right\rbrack}.$However, the packing is performed in a non-contiguous manner, such thattensor tiles 404 a-404 f do not cover each a contiguous area of matrix402. Thus, for example, tile tensor 404 a may include matrix 492elements 00, 02, 04, 30, 32, 34.FIGS. 4B and 4C show two exemplary 8-parallel kernel evaluations using a2×2 filter or kernel performed over matrix 402. The dashed rectangleillustrates the convolution operator of a [2,2] filter or kernel 403 onthe M[6,6] matrix 402. The tile tensors 404 a-404 f illustrate the sameoperation using tile tensor representation, using component-wisesummation of tiles. FIG. 4B illustrates the first convolution iteration,ignoring the empty slots. FIG. 4C illustrates a left circular rotationby 1.

Input tensor 402 is packed using tile tensors 404 a-404 f of shape[2,4], notated

${T_{M}\left\lbrack {\frac{6 \sim}{2},\frac{6 \sim}{4}} \right\rbrack}.$However, the packing is performed in a non-contiguous manner, such thattensor tiles 404 a-404 f do not cover each a contiguous area of matrix402. Thus, for example, tile tensor 404 a may include matrix 492elements 00, 02, 04, 30, 32, 34.

FIGS. 4B and 4C show two exemplary 8-parallel kernel evaluations using a2×2 filter or kernel performed over matrix 402. The dashed rectangleillustrates the convolution operator of a [2,2] filter or kernel 403 onthe M[6,6] matrix 402. The tile tensors 404 a-404 f illustrate the sameoperation using tile tensor representation, using component-wisesummation of tiles. FIG. 4B illustrates the first convolution iteration,ignoring the empty slots. FIG. 4C illustrates a left circular rotationby 1.

In some embodiments, interleaved tiling may be specified separately foreach dimension. For example, in

$\left\lbrack {\frac{5}{2},\frac{6 \sim}{4}} \right\rbrack$only the second dimension is interleaved. Also, although with basictiling it holds that

${e_{i} = \frac{n_{i}}{t_{i}}},$for interleaved tiling it is sometimes useful to have larger values fore_(i). In this case, this value can be explicitly stated using thenotation:

$\frac{n_{i} \sim e_{i}}{t_{i}}.$

FIGS. 5A-5C illustrate an interleaved tiling operation. In FIG. 5A,matrix 502 of size M[9,10] may be tiled using a tensor tile of size[2,4]. In FIG. 5A, a first tile tensor 504 a covers an area of matrix502 indicated by the dashed line rectangle, however, it does not coverthe area in a contiguous manner. Rather, first tile tensor 504 aincludes matrix elements indicated by solid line rectangles. FIGS. 5Band 5C show tile tensors 504 b and 504 c, respectively, which are spreadover matrix 502 using equal strides.

FIGS. 6A-6D also illustrate an interleaved tiling process of the presentdisclosure. As can be seen, input tensor 602 a is a region of matrix 602covering 3×5 elements. Input tensor 602 b may be packed into tiletensors 604, 606, 608. In FIG. 6B, input tensor 602 a is packed in aninterleaved manner into 15 tile tensors 604 a-604 o. the elements ofregion 602 a are packed into the first slot in tile tensors 604 a-604 o.In some embodiments, the same region may be packed in a duplicativemanner into all slots of the tile tensors (FIG. 6C). in otherembodiments, each tile tensor slot may pack a different region of matrix602 (FIG. 6D).

Convolution Using Tile Tensors

Background

As noted above, the present disclosure provides for a data packingscheme which permits efficient high-level tensor manipulationoperations, such as convolution, over the encrypted data, while reducingcomputational overhead.

Convolution may be described as a sliding window function applied to amatrix. FIG. 2A is a diagram of a basic convolutional neural networkarchitecture. FIG. 2B illustrates the sliding window (or filter, orkernel) function.

The input of a convolutional layer is often an images tensor, e.g.,image tensor 202 in FIGS. 2A and 2B, denoted I[w_(I), h_(I), c, b], anda filter tensor F[w_(F), h_(F), c, f], e.g., filter 204 in FIG. 2B withthe following shape parameters: width w_(I), W_(F), height h_(I), h_(F),and the number of image channels c (e.g., 3 for an RGB image). In somecases, convolution may be computed for a batch of b images and for ffilters. Informally, the convolution operator moves each filter in F asa sliding window over elements of I, starting at position (0,0) andusing strides of δ_(w) and δ_(h). When the filter fits entirely in theinput, the inner product is computed between the elements of the filterand the corresponding elements of I. For example, in FIG. 2B, inputmatrix 202 is a [5,5] matrix which may represent an image, wherein eachmatrix element corresponds to one pixel (e.g., a value between 0 and 255for grayscale images). The sliding window 204 may be referred to as thekernel, filter, or feature detector, and is a 3×3 filter. the filter 204slides over matrix 202 in equal strides of one pixel, multiply itsfilter values element-wise with the matrix 202 values, then sums them upto output convolved feature 206. To get the full convolution, this isperformed for each element by sliding the filter 204 over the wholematrix 202.

Let I[w_(I), h_(I), c, b] and F[w_(F), h_(F), c, f] be two input tensorsfor the convolution operator representing images and filters,respectively. The results of the operation O=conv2d(I, F) is the tensorO[w_(O), h_(O), f, b], where

${w_{O} = \frac{w_{I} - w_{F} + 1}{\delta_{w}}},{h_{O} = \frac{h_{I} - h_{F} + 1}{\delta_{h}}},$δ_(w) and δ_(h) are the strides and

$\begin{matrix}{{O\left( {i,j,m,n} \right)} = {{\sum}_{k = 0}^{w_{F} - 1}{\sum}_{l = 0}^{h_{F} - 1}{\sum}_{p = 0}^{c - 1}{I\left( {{i \cdot {+ k}},{j \cdot {+ l}},p,n} \right)}{{F\left( {k,l,p,m} \right)}.}}} & (5)\end{matrix}$

In the degenerated case where δ_(w)=δ_(h)=b=f=c=1, Eq(5) can besimplified to

$\begin{matrix}{{O\left( {i,j} \right)} = {{\sum}_{k = 0}^{w_{F} - 1}{\sum}_{l = 0}^{h_{F} - 1}{I\left( {{i + k},{j + l}} \right)}{{F\left( {k,l} \right)}.}}} & (6)\end{matrix}$

In FHE settings, it is sometimes useful to convert a convolutionoperation to a matrix-matrix multiplication by pre-processing the inputbefore encrypting it. One such method is image-to-column, which works asfollows for the case c=b=1. Given an image I[w_(I),h_(I)] and f filtersF[w_(F), h_(F), f], the operator I′, F′=im2col(I, F) computes a matrixI[w_(O)h_(O), w_(F)h_(F)], where each row holds the content of a validwindow location in I flattened to a row-vector, and F′[w_(F)h_(F), f]contains every filter of F flattened to a column-vector. Here, thetensor O′[w_(O)h_(O), f]=I′F′ is a flattened version of the convolutionresult O[w_(O), h_(O), f]=conv2d(I, F).

In some embodiments, the present disclosure provides for a variant I″,F″=im2col′(I, F) that computes I″[w_(O)h_(O)f, w_(F)h_(F)] byconsecutively replicating f times every row of I′, and F″[w_(O)h_(O)f,w_(F)h_(F)] by concatenating w_(O)h_(O) times the matrix F′^(T). Thetensor O″[w_(O)h_(O)f, 1]=sum(I″*F″, 2) contains the convolution resultO[w_(O), h_(O), f]. The advantage of this variant is that the output isfully flattened to a column vector, which is useful in situations whereflattening is costly (e.g., in FHE). The drawback of this variant methodis that it is impossible to perform two consecutive convolutionoperators without costly pre-processing in between.

In some embodiments, the present disclosure provides for a novel methodto compute convolution over FHE data. In some embodiments, the presentmethod provides for greater efficiency when the input is a large image,and allows for efficient computation of consecutive convolution layersin a FHE only system.

Convolution with Interleaved Tile Tensors

In some embodiments, the present disclosure provides for usinginterleaved tile tensors, as detailed above, to efficiently computeconvolution over FHE data.

With reference back to FIGS. 4A-4C, matrix M[6,6] 402 is packed using 6different tile tensors 404 a-404 f denoted

${T_{M}\left\lbrack {\frac{6 \sim}{2},\frac{6 \sim}{4}} \right\rbrack}.$Here, the tile shape is [2,4] and the external tensor shape, i.e., thetotal number and arrangement of tile tensors used to pack matrix 402 is[3,2]. Every tile contains a 2×4 sub-matrix, but instead of beingcontiguous, it comprises a set of matrix elements spaced evenly in thematrix. For example, as can be seen in FIG. 4B, the elements00,01,10,11,20,21 are all mapped to the same slot in the different tiletensors 404 a-404 f. In other words, the matrix 402 region indicated bya dashed line rectangle 403 is all mapped to the upper left slot in eachof tile tensors 404 a-404 f.

The interleaved packing allows for a more efficient implementation ofEquation 6 above with respect to runtime and storage. Intuitively, SIMDwas used to compute multiple elements of the output in a singleoperation. The filter is packed simply as

${T_{F}\left\lbrack {w_{F},h_{F},{\frac{*}{t_{1}}{,\frac{*}{t_{2}}}}} \right\rbrack},$i.e., it has w_(F)h_(F) tiles, each containing one value of the filterin all slots. This allows multiplying each image tile with each value ofthe filter.

For example, FIG. 4B shows a computation of the convolution output whena [2,2] filter 403, indicated by dashed line rectangle, is placed at thetop left position of matrix 402. The SIMD nature of the computationcomputes the output in other regions as well. The result is a singletile, where each slot contains the convolution result of thecorresponding region, such that this tile is packed in the sameinterleaved packing scheme as the input tiles.

A more complicated example is given in FIG. 4C. Here the filter 403 isplaced one pixel to the right. As a result, the filter 403 needs to bemultiplied by elements that appear in different regions, i.e. they aremapped to slots of different indices. In this case, the tiles need to berotated appropriately. For example, placing the filter with its upperleft corner on pixel (01), the convolution is computed using the (0,0)slot of tiles 404 d and 404 e, and slot (0,1) of tiles 404 a and 404 b.accordingly, tiles 404 a and 404 b are therefore rotated to move therequired value to slot (0,0) as well.

The total cost of convolution when using this packing is summarized inthe following lemma:

-   -   Let s be the number of slots in a ciphertext. Then, given an        input image I[w_(I), h_(I)] and a filter F[w_(F), h_(F)],        packing I as

$T_{I}\left\lbrack {\frac{w_{I} \sim}{t_{1}},\frac{h_{I} \sim}{t_{2}}} \right\rbrack$

-   -    and the filter as

${T_{F}\left\lbrack {w_{F},h_{F},{\frac{*}{t_{1}}{,\frac{*}{t_{2}}}}} \right\rbrack},$

-   -    convolution can be computed using: O(w_(I)h_(I)w_(F)h_(F)/s)        multiplications, and

$O\left( {{w_{F}\frac{w_{I}}{t_{1}}} + {h_{F}\frac{h_{I}}{t_{2}}} + {w_{F}h_{F}}} \right)$

-   -    rotations. The input is encoded in O(w_(I)h_(I)/s) ciphertexts.    -   Proof: Multiplications. To compute the convolution, each of the        w_(I)h_(I) elements of the input tensor must be multiplied with        each of the w_(F)h_(F) elements of the filter, excluding edge        cases that do not change the asymptotic behavior. Since each        multiplication multiplies s slots, only        O(w_(I)h_(I)w_(F)h_(F)/s) multiplications are needed.    -   Rotations. Recall the output is of size        (w_(I)−w_(F)+1)(h_(I)−h_(F)+1). The k-th slot of different        ciphertexts is mapped to elements of I with indexes

${k\frac{w_{I}}{t_{1}}} \leq x_{o} < {\left( {k + 1} \right)\frac{w_{I}}{t_{1}}{and}k\frac{h_{I}}{t_{2}}} \leq y_{o} < {\left( {k + 1} \right){\frac{h_{I}}{t_{2}}.}}$

-   -    It is therefore enough to analyze the cost of computing the        convolution for

${0 \leq x_{o} < {\frac{w_{O}}{t_{1}}{and}0} \leq y_{o} < \frac{h_{O}}{t_{2}}},$

-   -    since computing the other elements of the output has no cost        due to the SIMD feature. It follows that a rotation is needed        when

${x_{o} + i} \geq {{\frac{w_{I}}{t_{1}}{or}y_{o}} + j} \geq {\frac{h_{I}}{t_{2}}.}$

-   -    This totals to

${O\left( {{w_{F}\frac{w_{I}}{t_{1}}} + {h_{F}\frac{h_{I}}{t_{2}}} + {w_{F}h_{F}}} \right)}.$

-   -   Storage. Because O(s) slots of each ciphertext are used, the        input can be encoded in O(w_(I)h_(I)/s) ciphertexts.

The output of the convolution is a tile tensor

${T_{O}\left\lbrack {\frac{w_{O} \sim ?}{t_{1}},\frac{h_{O} \sim ?}{t_{2}}} \right\rbrack}.$The unknown values are introduced by filter positions that extend beyondthe image, as shown in FIG. 4C. Note further that the external sizes

$e_{1} = {{\frac{w_{I}}{t_{1}}{and}e_{2}} = \frac{h_{I}}{t_{2}}}$of the tile tensor T_(I) remain the same in T_(O), and they may belarger than those actually required to hold the tensor O[w_(O), h_(O)].Hence, a more accurate depiction of T_(O)'s shape is

${T_{O}\left\lbrack {\frac{w_{O} \sim {e_{1}?}}{t_{1}},\frac{h_{O} \sim {e_{2}?}}{t_{2}}} \right\rbrack},$however, this may be ignored in practice.Handling Strides, Batching and Multiple Channels and Filters

In some embodiments, the present convolution algorithm may be extendedto handle multiple channels (e.g., images having multiple colorchannels, such as RGB images having red, green, and blue channels),multiple filters, and batching (e.g., processing multiple imagestogether).

Let a received input data be a tensor of images I[w_(I), h_(I), c, b],where c is the number of channels and b is the batch size. The input isthen packed as

${T_{I}\left\lbrack {\frac{w_{I} \sim}{t_{1}},\frac{h_{I} \sim}{t_{2}},\frac{c}{t_{3}},\frac{b}{t_{4}},\frac{*}{t_{5}}} \right\rbrack},$together with filters F[w_(F), h_(F), c, f], where f is the number offilters, as

${T_{F}\left\lbrack {w_{F},h_{F},{\frac{*}{t_{1}}{,{\frac{*}{t_{2}}{,\frac{c}{t_{3}},{\frac{*}{t_{4}}{,\frac{f}{t_{5}}}}}}}}} \right\rbrack},$where t_(i)∈

and Πt_(i)=s.

The convolution is computed similarly to the description above,multiplying tiles of T_(I) with the appropriate tiles of T_(F). Theresult is a tile tensor of shape

${T_{O}\left\lbrack {\frac{w_{O^{\sim}}?}{t_{1}},\frac{h_{O^{\sim}}?}{t_{2}},\frac{c}{t_{3}},\frac{b}{t_{4}},\frac{f}{t_{5}}} \right\rbrack}.$Summing over the channel (i.e., third) dimension using

$O\left( {\frac{w_{O}}{t_{1}}\frac{h_{O}}{t_{2}}\frac{b}{t_{4}}\frac{f}{t_{5}}\log t_{3}} \right)$rotations, to obtain

${T_{O}\left\lbrack {\frac{w_{O^{\sim}}?}{t_{1}},\frac{h_{O^{\sim}}?}{t_{2}},\frac{1?}{t_{3}},\frac{b}{t_{4}},\frac{f}{t_{5}}} \right\rbrack}.$

For bigger strides, >1 (resp. >1), it is required that either

$t_{1} = {{1\left( {{{resp}.t_{2}} = 1} \right){or}\frac{h_{I}}{t_{2}}{mod}} = {0{\left( {{{{resp}.\frac{w_{I}}{t_{1}}}{mod}} = 0} \right).}}}$Then, the implementation trivially skips ciphertexts in every row andciphertexts in every column.A Sequence of Convolutions

In some embodiments, the present disclosure provides for implementing asequence of multiple convolution layers. This is something that iscommon in neural networks. One of the advantages of the tile tensormethod is that the output of one convolution layer can be easilyadjusted to be the input of the next convolution layer.

In some embodiments, an input may be a batch tensor I[w_(I), h_(I), c,b] and a sequence of convolution layers with the l'th layer having afilter tensor F^(l)[w_(F) ^(l), h_(F) ^(l), c^(l), f^(l)]. For the firstlayer, c¹=c, and, for subsequent layers l>1, c^(l)=f^(l−1). As before,the input tensor may be packed as

${T_{I}\left\lbrack {\frac{w_{I^{\sim}}}{t_{1}},\frac{h_{I^{\sim}}}{t_{2}},\frac{c}{t_{3}},\frac{b}{t_{4}},\frac{*}{t_{5}}} \right\rbrack}.$For odd layers, l=2l+1, the filter tensor may be packed as

${T_{F}^{l}\left\lbrack {w_{F}^{l},h_{F}^{l},{\frac{*}{t_{1}}{,{\frac{*}{t_{2}}{,\frac{c}{t_{3}},{\frac{*}{t_{4}}{,\frac{f^{l}}{t_{5}}}}}}}}} \right\rbrack}.$The output is then

${T_{O}\left\lbrack {\frac{w_{O^{\sim}}^{l}?}{t_{1}},\frac{h_{O^{\sim}}^{l}?}{t_{2}},\frac{1?}{t_{3}},\frac{b}{t_{4}},\frac{f^{l}}{t_{5}}} \right\rbrack}.$For even layers, l=2l, the filters may be packed as:

${T_{F}^{l}\left\lbrack {w_{F}^{l},h_{F}^{l},{\frac{*}{t_{1}}{,{\frac{*}{t_{2}}{,\frac{f}{t_{3}},{\frac{*}{t_{4}}{,\frac{c}{t_{5}}}}}}}}} \right\rbrack}.$

As can be seen, the shapes of the layer outputs do not match the shapesof the inputs of the subsequent layers. Accordingly, in someembodiments, the present disclosure provides for adjusting an output ofan odd layer to be suitable for the next even layer. For this, in someembodiments, the present disclosure provides for clearing the unknownsby multiplying with a mask, and then replicating the channel dimension.The result is a tile tensor of this shape:

${T_{O}\left\lbrack {\frac{w_{O^{\sim}}^{l}?}{t_{1}},\frac{h_{O^{\sim}}^{l}?}{t_{2}},{\frac{*}{t_{3}}{,\frac{b}{t_{4}},\frac{f^{l}}{t_{5}}}}} \right\rbrack},$which matches the input format of the next layer, because f^(l)=c^(l+1).To make an output of an even layer suitable for the next odd layer, thepresent disclosure similarly cleans and replicates along the filterdimension.

It may be noted that changing the order of the dimensions leads to asmall improvement. The improvement comes because summing over the firstdimension ends up with a replication over this dimension. Therefore,setting the channel dimension as the first dimension saves thereplication step when preparing the input to an even layer. In someembodiments, cleaning may be skipped as well, because the unknown valuesalong the image width and height dimensions do not affect the result.Alternatively, the filter dimension can be set as first and then thereplication step can be skipped when preparing the input for an oddlayer.

Naive Convolution Methods The above method reduces to a simple methodknown by various names such as SIMD packing when t₁=t₂=t₃=t₅=1. In thiscase, every element in the tensors for the images and filters is storedin a separate ciphertext, and the slots are only used for batching. Insome embodiments, the reduction to matrix multiplication as describedabove may also be used, however, it may only be applicable in the caseof neural networks with one convolutional layer.

Experimental Results

The present inventors conducted experiments using two benchmark neuralnetwork models:

-   -   CryptoNets (see, Ran Gilad-Bachrach, et al. Cryptonets: Applying        neural networks to encrypted data with high throughput and        accuracy. In International Conference on Machine Learning, pages        201-210, 2016.). This network has a convolutional layer followed        by two fully connected layers.    -   AlexNet AlexNet (see, Alex Krizhevsky et al. Imagenet        classification with deep convolutional neural networks. Neural        Information Processing Systems, 25, 01 2012.)

The models were trained on the MNIST (see, Yann LeCun, et al. The MNISTdatabase of handwritten digits. 10:34, 1998) and COVIDx CT-2A (see,Hayden Gunraj, et al. Covidnet ct-2: Enhanced deep neural networks fordetection of covid-19 from chest CT images through bigger, more diverselearning. arXiv preprint arXiv:2101.07433, 2021) data-sets,respectively.

The results are reported herein with respect to performing modelinference using these model weights in encrypted and unencrypted forms.AlexNet was sued to demonstrate the power of the present disclosure, andCryptoNets was used to demonstrate the effect of different packing onthe computation performance and memory.

The experiments were run using a computer running an Intel Xeon CPUE5-2699 v4 @ 2.20 GHz, with 44 cores (88 threads) and 750 GB memory. Theexperiments used the CKKS SEAL implementation targeting 128 bitssecurity, and all the reported results are the average of at least 10runs.

CryptoNets Results

The CryptoNets model was implemented using tile tensors of shape

$\left\lbrack {\frac{n_{1}}{t_{1}},\frac{n_{2}}{t_{2}},\frac{b}{t_{3}}} \right\rbrack,$where b is the batch size. In practice, reported are only the resultsfor the case t₃=b that minimizes the overall latency by filling all theciphertext slots (8192 in this case). For the convolution layer, thenaïve SIMD method was used, when b equals the number of plaintext slotsand t₁=t₂=1. Otherwise, the present variant of the im2col operator wasused (see above). These methods work better than the present novelconvolution operator when the images are small and the network has oneconvolutional layer.

FIG. 7 is an illustration of the present CryptoNets implementation usingtile tensors. For simplicity, only the component-wise square activationlayers are indicated by the sq( ) function because they maintain thetile tensor shape. The equations on the right represent the underlyingtensor operations. The input tensors are I, F, B_(c), W₁, B₁, W₂, B₂,where I′, F′=im2col′(I, F). FIG. 7 shows the tile tensor flow in thepresent implementation. Here, the inputs I and F are the image andfilter matrices, respectively, and I′, F′=im2col′(I, F). In addition,B_(c) is the trained bias of the convolution layer and W₁, W₂, B₁, B₂are the trained weights and bias info of the FC layers.

Table 1 reports the latency and memory usage for performing a modelinference with different tile shapes when t₃=b=1. For brevity, t₁ isonly considered to be at the extreme points (e.g., t₁=1,8192) or t₁values that led to optimal solution, and some additional samples. Theoptimal latency and memory usage are achieved for t₁=32, which allowspacking the tensors I, F, W₁ using the minimal number of tiles.

Table 2 reports the latency, amortized latency, and memory usage forperforming a model inference with different t₃=b values. For every suchvalue, reported are only the t₁, t₂ values that led to the optimalsolutions. Unlike the case where b=1, here every choice of t₃ leads to adifferent trade-off between the performance measures. For example, whenincreasing t₃, the latency and memory consumption increase, but theper-sample amortized latency decreases. The encryption and decryptiontime also increase with t₃, except for the case t₃=8192, where the naiveSIMD convolution operator is used.

TABLE 1 Running a model inference with different tile shapes [t₁, t₂,t₃] when t₃ = b = 1. Latency Enc + Dec Memory t₁ t₂ t₃ (sec.) (sec.)(GB) 1 8192 1 0.86 0.04 1.58 8 1024 1 0.56 0.04 0.76 32 256 1 0.56 0.040.73 64 128 1 0.57 0.04 0.77 128 64 1 0.61 0.04 0.94 256 32 1 0.68 0.051.37 1024 8 1 1.93 0.14 3.17 8192 1 1 11.10 0.80 14.81 The reportedvalues are the inference latency, the encryption and decryption time,and the memory usage peak.

TABLE 2 Amortized Latency Latency Enc + Dec Memory t₁ t₂ t₃ (sec.)(sec.) (sec.) (GB) 32 256 1 0.56 0.56 0.04 0.73 16 128 4 0.56 0.14 0.051.20 8 64 16 0.6 0.037 0.10 2.49 4 32 64 0.95 0.015 0.24 6.62 1 32 2561.94 0.008 0.70 16.38 1 8 1024 5.6 0.0055 2.68 61.45 1 2 4096 21.570.0053 12.55 242.46 1 1 8192 41.32 0.005 1.29 354.47 Running a modelinference with different tile shapes [t₁, t₂, t₃], reporting only theoptimal t₁ and t₂ choices for a range of different t₃ = b values. Thereported values are: the inference latency, the amortized latency(latency/b), the encryption and decryption time, and the memory usagepeak.AlexNet Benchmark

For the AlexNet benchmark, a variant of AlexNet network is used, whichincludes 5 convolution layers, 3 fully connected layers, 7 ReLUactivations, 3 Batch Normalization layers, and 3 Max Pooling layers. ACKKS-compliant variant of AlexNet was created by replacing the ReLU andMax Pooling components with a scaled square activation and AveragePooling correspondingly along with some additional changes. It wastrained and tested on the COVIDx CT-2A dataset. The COVIDx CT-2A datasetis an open access benchmark of CT images dataset that contains threeclasses of chest CT images: Normal, Pneumonia, and COVID-19 cases. Theexperiments used a subset of 10,000 images per class for training, 1,000images per class for validation, and 201 images in total for test with67 random samples from each class. The images were resized to 224*224*3to fit the input size expected by AlexNet.

The biases were packed in 5-dimensional tile tensors with compatibleshapes, allowing us to add them to the convolution outputs. The fullyconnected layers were handled using the matrix-matrix multiplicationtechnique (see above). The input to these layers arrives from theconvolutional layers as a 5-dimensional tile tensor,

$\left\lbrack {\frac{*}{t_{1}}{,\frac{1 \sim}{t_{2}},\frac{1 \sim}{t_{3}},\frac{256}{t_{4}},\frac{b}{t_{5}}}} \right\rbrack.$Therefore, the first fully connected layer is packed in 5 dimensions aswell:

$\left\lbrack {\frac{4096}{t_{1}},\frac{1 \sim}{t_{2}},\frac{1 \sim}{t_{3}},\frac{256}{t_{4}},\frac{b}{t_{5}}} \right\rbrack.$Its output,

$\left\lbrack {\frac{4096}{t_{1}},\frac{1 \sim}{t_{2}},\frac{1 \sim}{t_{3}},\frac{1?}{t_{4}},\frac{b}{t_{5}}} \right\rbrack,$is replicated along dimensions 2 through 4, then flattened using theflatten operator to

$\left\lbrack {\frac{4096}{t_{1}},{\frac{*}{t_{2}t_{3}t_{4}}{,\frac{b}{t_{5}}}}} \right\rbrack,$from which it may continue normally.

The accuracy of running regular AlexNet and the HE-friendly AlexNet weremeasured using PyTorch (PyTorch library 1.5.1, see https://pytorch.org)over a plaintext test-set. The results were 0.861 and 0.806,respectively. No additional accuracy degradation was observed whenrunning the HE-friendly AlexNet using the present framework overencrypted data.

Table 3 reports the time and memory consumption for the latterexperiment using 4 configurations on a set of 30 representative samples.The configurations involve unencrypted model weights (PT) and encryptedmodel weights (CT) optimized for low latency (Latency) or highthroughput (TP). For these configurations, the inference results werealso compared with the inference results of running HE-Friendly AlexNeton PyTorch over the plaintext test-set by calculating the Root MeanSquare Error (RMSE). These were always less than 4e−3.

TABLE 3 AlexNet executed in the present framework with differentconfigurations. Latency Amortized Enc + Dec Memory Configuration (sec.)Latency (sec.) (sec.) (GB) PT-Latency 181.9 181.9 5.3 123.8 PT-TP 720.890.1 5.4 568.1 CT-Latency 358.1 358.1 5.4 223.4 CT-TP 1130.4 282.6 5.6688.8Optimizer Accuracy

The optimizer module 110 (see FIGS. 1A, 1B) and its simulator 110 bestimate the time and memory usage for a given configuration option on asingle CPU thread. For that, it relies on pre-benchmarked measures ofthe different FHE operations. To assess the accuracy of theseestimations, the following experiment on HE-friendly AlexNet wereperformed using encrypted model. The four configuration options wereselected which achieved the lowest estimated latency when using localsearch and compared the inference time and the encryption time of theinput and the model between the simulation output and an actual run overencrypted data. Table 4 summarizes the results. It may be observed thatthe simulator provides relatively accurate time estimations for all fourconfigurations. The average estimated time deviation is −15.8%, −11.9%,and −7.2% for inference, model encryption, and batch input encryption,respectively. It is noted that the simulated storage matches themeasured storage for all configurations, thus this data is not includedin Table 4.

TABLE 4 Inference Model encryption. Input encryption. Configuration(sec.) (sec.) (sec.) Conf₁ 4232 (−11%)   1509 (−11.5%) 162 (−6.8%) Conf₂4758 (−13.9%) 1493 (−12.1%) 164 (−7.9%) Conf₃ 4927 (−18.1%) 1680(−11.5%) 177 (−6.8%) Conf₄ 4798 (−20%)   1668 (−12.3%) 178 (−7.3%)Simulated time estimations for the configurations (Conf_(i))_(i=1..4)formatted as [tile shape − convolution mode]: [16, 8, 8, 16, 1]-CWHFB,[8, 8, 8, 32, 1]-CWHFB, [16, 8, 8, 16, 1]-FWHCB, [32, 8, 8, 8, 1]-FWHCB,respectively. The acronyms CWHFB and FWHCB indicate the order ofdimensions in the tile tensor. The deviation of the estimated times fromthe real times are reported in brackets.Additional Comparisons

Tile tensors capture as a special case the simple method where eachelement of the input matrices is placed in a separate ciphertext. Table2 reports the results for this method in the last row.

Two more special cases of matrix-vector multiplication algorithms aredescribed in [Crockett 2020] (see, Eric Crockett. A low-depthhomomorphic circuit for logistic regression model training. CryptologyePrint Archive, Report 2020/1483, 2020.) These are equivalent to Eq(1)and Eq(2) above. In addition, [Crockett 2020] shows an extension tomatrix-matrix multiplication by extracting columns from the secondmatrix and applying matrix-vector multiplication with each. Thisextraction of columns requires multiplication by mask and increases themultiplication depth. With the present tile tensors method, a naturalextension to matrix-matrix multiplication may be obtained that doesn'trequire increasing the multiplication depth.

A different family of techniques are based on diagonalization. The basicmethod for matrix-vector multiplication is described in [Halevi 2014](see, Shai Halevi and Victor Shoup. Algorithms in helib. In Juan A.Garay and Rosario Gennaro, editors, Advances in Cryptology—CRYPTO2014-34th Annual Cryptology Conference, Santa Barbara, CA, USA, Aug.17-21, 2014, Proceedings, Part I, volume 8616 of Lecture Notes inComputer Science, pages 554-571. Springer, 2014.doi:10.1007/978-3-662-44371-2_31.) For a ciphertext with n slots, an n×nmatrix is pre-processed to form a new matrix, where each row is adiagonal of the original matrix. Then, multiplication with a vector canbe done using n rotations, multiplications, and additions. Theperformance of the present depends on the tile shape. For example, forsquare tiles of a shape approximating [√{square root over (n)}, √{squareroot over (n)}], the matrix-vector multiplication costs nmultiplications and √{square root over (n)} log √{square root over (n)}rotations. (The matrix breaks down to n tiles in this case; each needsto be multiplied with one vector tile. The summation reduces the shapeof the external tensor to [√{square root over (n)}, 1], and each of theremaining tiles is summed over using log √{square root over (n)}rotation).

Some improvements to diagonalization techniques have been presented,which reduce the number of required rotations to O(√{square root over(n)}) under some conditions, and by exploiting specific properties ofthe HE schemes of HElib [Halevi 2014]. The present methods make nospecial assumptions. Exploiting such properties and combining them withthe tile tensor data structure is reserved for future work.

In [Ciaoqian 2018] (see, Xiaoqian Jiang, et al. Secure outsourced matrixcomputation and application to neural networks. In Proceedings of the2018 ACM SIGSAC Conference on Computer and Communications Security, CCS'18, page 1209-1222, New York, NY, USA, 2018. Association for ComputingMachinery.), a matrix-matrix multiplication method based ondiagonalization is described. They reduce the number of rotations toO(n) instead of O(n₂) for multiplying with n vectors. However, thiscomes at the cost of increasing the multiplication depth by twomultiplications with plaintexts. This is a significant disadvantage innon-client-aided FHE, since the performance of a circuit is generallyquadratic in its depth, and from practical considerations the depth issometimes bounded.

Convolution

A convolution layer is a basic building block in NN. Previous workpresented optimizations for small input: GAZELLE (see, Chiraag Juvekar,et al. GAZELLE: A low latency framework for secure neural networkinference. In 27th USENIX Security Symposium (USENIX Security 18), pages1651-1669, Baltimore, MD, August 2018. USENIX Association.) considered a28×28 grey scale images. GALA (see, Qiao Zhang, et al. Gala: Greedycomputation for linear algebra in privacy-preserved neural networks.arXiv preprint arXiv:2105.01827, 2021.) considered 16×16 images. HEAR(see, Miran Kim, et al. HEAR: human action recognition via neuralnetworks on homomorphically encrypted data. CoRR, abs/2104.09164, 2021.)considered 3D tensor input of size 32×15×2.

In contrast, the present inventors considered 224×224 RGB images. Usingthe methods of Gazelle, GALA, and/or HEAR for such large inputs is lessefficient, because they pack c_(n) channels of the input in a singleciphertext. Then, they act on all c_(n) channels taking advantage of theSIMD feature. For example, GALA and GAZELLE require a total of

$O\left( \frac{f + {{cw}_{I}h_{I}}}{c_{n}} \right)$rotation and multiplication operations, where f, c, w_(I), h_(I) aredefined above. In HEAR, a sequence of convolutions is considered. Then apre-processing step between two convolution steps is needed. Computingthe pre-processing step and the convolution takes

$O\left( {w_{F}h_{F}\frac{w_{I}}{t_{1}}\frac{h_{I}}{t_{2}}\frac{c^{2}}{t_{3}}f} \right)$rotation and multiplication operations. For images of size244×244=50,176 at most one channel can fit in a ciphertext that has65,536 slots, i.e. c_(n)=1. Using ciphertexts with fewer slots or biggerimages results in performance degradation since the data of a singlechannel is spread among several ciphertexts. Previous work did notexplain how to extend to support this case efficiently. Trivially, moreslots can be emulated using several ciphertexts. This adds to therunning time a factor proportional to the image size, i.e.O(w_(I)h_(I)).

In the present convolution method, the number of rotations for largerimages increases by a factor of

$O\left( {{w_{F}\frac{w_{I}}{t_{1}}} + {h_{F}\frac{h_{I}}{t_{2}}} + {w_{F}h_{F}}} \right)$and the number of multiplications by

${O\left( {\frac{w_{I}}{t_{1}}\frac{h_{I}}{t_{2}}} \right)},$which is better than previous works for large images. For multiplechannels, filters and samples in a batch, the present method run timeincreases by a factor of

${O\left( {\frac{c}{t_{3}}\frac{b}{t_{4}}\frac{f}{t_{5}}} \right)},$and additional

$O\left( {\frac{w_{O}}{t_{1}}\frac{h_{O}}{t_{2}}\frac{b}{t_{4}}\frac{f}{t_{5}}\log t_{3}} \right)$rotations required for summing over channels inside the tile. Bychoosing values for the tile shape t_(i), it may be optimized for theparticular sizes of a given computation.Sequence of Convolution Layers

In GAZELLE and GALA, optimizations were made for a single convolutionlayer. While this is important, deep networks have long sequences ofconvolution networks of different sizes and different numbers offilters. For example, AlexNet has five consecutive layers of convolutionof different sizes.

To support more layers, previous works assumed a non FHE step, such asgarbled circuits or another MPC protocol, after each layer (in aclient-aided approach). The non-FHE step performs the activationfunction and puts the input for the next layer in the correct packing.Converting the packing using an FHE-only system is expensive. In HEAR,an all-FHE solution was considered. However, they required apre-processing step that needs O(w_(I)h_(I)cb) multiplications and

$O\left( {w_{I}h_{I}{cb}\frac{{t_{1}t_{2}} - 1}{t_{2}t_{2}}} \right)$rotations.

In contrast, the present packing method requires a pre-processing stepbefore even layers only. In that case, it requires

$O\left( {\frac{w_{I}}{t_{1}}\frac{h_{I}}{t_{2}}\frac{c}{t_{3}}\frac{b}{t_{4}}\log t_{5}} \right)$rotations; here, w_(I), h_(I), c and b refer to the image dimensions,the number of channels and the batch size in the input to the layer.Neural Network Inference

The present approach is comrade with other end-to-end neural networkinference solutions, e.g., nGraph-HE2 (see, Fabian Boemer, et al.NGraph-HE2: A High-Throughput Framework for Neural Network Inference onEncrypted Data. In Proceedings of the 7^(th) ACM Workshop on EncryptedComputing & Applied Homomorphic Cryptography, WAHC'19, pages 45-56, NewYork, NY, USA, 2019. Association for Computing Machinery.) and TenSEAL(see, Ayoub Benaissa, et al. TenSEAL: A Library for Encrypted TensorOperations Using Homomorphic Encryption. arXiv, 2021.)

Table 5 reports the comparison results. TenSEAL uses diagonalizationtechniques for matrix-multiplication and im2col for convolution,assuming a single image as input. Moreover, TenSEAL assume unencryptedmodel weights. Hence, TenSEAL was comrade to the present framework whenoptimized for batch size of one, for unencrypted model weights (PT) andfor completeness also show results for encrypted model weights (CT).nGraph-HE2 also focuses on unencrypted models. It uses SIMD packing,which is a special case of the present framework when optimized for thelargest possible batch size.

TABLE 5 Comparison of the CryptoNets benchmark with the present tiletensor framework and other NN compilers that are freely availableonline. Latency Amortized Framework (sec.) Latency (sec.) TenSeal (b= 1) 3.55 3.55 Present Method-PT (b = 1) 0.48 0.48 Present Method-CT (b= 1) 0.56 0.56 nGraph-HE2 (b = 8192) 11.93 0.00146 Present Method-TP (b= 8192) 13.52 0.00165 Present Method-CT (b = 8192) 41.32 0.00504 b = 1and b = 8192 were set for the top and bottom lines, respectively.

The results highlight the efficiency and versatility of the presentframework. Targeting optimal latency, the present framework provides atleast seven times speed-up over nGraph-HE2 and TenSEAL. Moreover, it canadapt to variable batch sizes. When targeting optimal throughput,nGraph-HE2 was slightly faster than the present framework. This can beexplained by the fact that the present library currently focuses onoptimizing the packing scheme, which in this case are identical to theone used by nGraph-HE2. Hence, the two libraries perform the exact sameset of homomorphic operations, but nGraph-HE2 also providesoptimizations for pipelining the underlying FHE instructions (e.g., bylazy rescaling). It may be stressed that the power of using differentpacking schemes is more noticeable for large networks that involve asequence of operations and is often not reflected in small networks suchas CryptoNets.

An additional framework that is not included in the above comparisonexperiments is the CHET compiler (see, Roshan Dathathri, et al. Chet: Anoptimizing compiler for fully-homomorphic neural-network inferencing. InProceedings of the 40th ACM SIGPLAN Conference on Programming LanguageDesign and Implementation, PLDI 2019, page 142-156, New York, NY, USA,2019. Association for Computing Machinery.), which performs inference ofencrypted data in a non-encrypted network. They report 2.5 secondslatency on a similarly sized, though less accurate, MNIST neural networkclassifier using 16 threads. They use a similar approach of an abstractdata structure, CipherTensor, combined with automatic optimizations. Thepresent believe CipherTensors are less flexible than tile tensors. Theyinclude a small fixed set of implemented layouts, each with its ownkernel of algorithms, whereas tile tensors offer a wider variety ofoptions with a single set of generalized algorithms. Further, it wasn'tdemonstrated that CipherTensors offer an easy method to trade latencyfor throughput and control memory consumption, as is possible in tiletensors by controlling the batch dimension. Finally, CipherTensorsrequire replication of the input data using rotations, whereas some ofthese replications can be avoided using tile tensors.

The EVA compiler (see, Roshan Dathathri, et al. Eva: An encrypted vectorarithmetic language and compiler for efficient homomorphic computation.In Proceedings of the 41st ACM SIGPLAN Conference on ProgrammingLanguage Design and Implementation, PLDI 2020, page 546-561, New York,NY, USA, 2020. Association for Computing Machinery.) is built on top ofCHET. They report an improved performance of 0.6 seconds on the samenetwork using 56 threads and various optimizations unrelated to packing;these optimizations are outside the scope of this paper. The presentbest result of 0.48 seconds was achieved for the more accurateCryptoNets architecture. The present inventors believe even betterresults can be obtained by combining the present packing optimizationswith EVA's optimizations (e.g., eliminating rescale operations to reducethe overall prime chain length).

The LoLa network (see, Alon Brutzkus, et al. Low latency privacypreserving inference. In Kamalika Chaudhuri and Ruslan Salakhutdinov,editors, Proceedings of the 36th International Conference on MachineLearning, volume 97 of Proceedings of Machine Learning Research, pages812-821, Long Beach, California, USA, 09-15 Jun. 2019. PMLR.) alsoreports results for the CryptoNets architecture. They achieve a latencyof 2.2 seconds using 8 threads. The LoLa network uses 150ciphertext-ciphertext multiplications, 279 rotations, and 399 additionsfor a single prediction. (the present deduced these numbers from LoLa'sdetailed description.) The present approach requires 32 multiplications,89 rotations, and 113 additions. These differences roughly explain theobserved latency results.

The present invention may be a system, a method, and/or a computerprogram product. The computer program product may include a computerreadable storage medium (or media) having computer readable programinstructions thereon for causing a processor to carry out aspects of thepresent invention.

The computer readable storage medium can be a tangible device that canretain and store instructions for use by an instruction executiondevice. The computer readable storage medium may be, for example, but isnot limited to, an electronic storage device, a magnetic storage device,an optical storage device, an electromagnetic storage device, asemiconductor storage device, or any suitable combination of theforegoing. A non-exhaustive list of more specific examples of thecomputer readable storage medium includes the following: a portablecomputer diskette, a hard disk, a random access memory (RAM), aread-only memory (ROM), an erasable programmable read-only memory (EPROMor Flash memory), a static random access memory (SRAM), a portablecompact disc read-only memory (CD-ROM), a digital versatile disk (DVD),a memory stick, a floppy disk, a mechanically encoded device havinginstructions recorded thereon, and any suitable combination of theforegoing. A computer readable storage medium, as used herein, is not tobe construed as being transitory signals per se, such as radio waves orother freely propagating electromagnetic waves, electromagnetic wavespropagating through a waveguide or other transmission media (e.g., lightpulses passing through a fiber-optic cable), or electrical signalstransmitted through a wire. Rather, the computer readable storage mediumis a non-transient (i.e., not-volatile) medium.

Computer readable program instructions described herein can bedownloaded to respective computing/processing devices from a computerreadable storage medium or to an external computer or external storagedevice via a network, for example, the Internet, a local area network, awide area network and/or a wireless network. The network may comprisecopper transmission cables, optical transmission fibers, wirelesstransmission, routers, firewalls, switches, gateway computers and/oredge servers. A network adapter card or network interface in eachcomputing/processing device receives computer readable programinstructions from the network and forwards the computer readable programinstructions for storage in a computer readable storage medium withinthe respective computing/processing device.

Computer readable program instructions for carrying out operations ofthe present invention may be assembler instructions,instruction-set-architecture (ISA) instructions, machine instructions,machine dependent instructions, microcode, firmware instructions,state-setting data, or either source code or object code written in anycombination of one or more programming languages, including anobject-oriented programming language such as Java, Smalltalk, C++ or thelike, and conventional procedural programming languages, such as the “C”programming language or similar programming languages. The computerreadable program instructions may execute entirely on the user'scomputer, partly on the user's computer, as a stand-alone softwarepackage, partly on the user's computer and partly on a remote computeror entirely on the remote computer or server. In the latter scenario,the remote computer may be connected to the user's computer through anytype of network, including a local area network (LAN) or a wide areanetwork (WAN), or the connection may be made to an external computer(for example, through the Internet using an Internet Service Provider).In some embodiments, electronic circuitry including, for example,programmable logic circuitry, a field-programmable gate array (FPGA), ora programmable logic array (PLA) may execute the computer readableprogram instructions by utilizing state information of the computerreadable program instructions to personalize the electronic circuitry,in order to perform aspects of the present invention. In someembodiments, electronic circuitry including, for example, anapplication-specific integrated circuit (ASIC), may be incorporate thecomputer readable program instructions already at time of fabrication,such that the ASIC is configured to execute these instructions withoutprogramming.

Aspects of the present invention are described herein with reference toflowchart illustrations and/or block diagrams of methods, apparatus(systems), and computer program products according to embodiments of theinvention. It will be understood that each block of the flowchartillustrations and/or block diagrams, and combinations of blocks in theflowchart illustrations and/or block diagrams, can be implemented bycomputer readable program instructions.

These computer readable program instructions may be provided to aprocessor of a general-purpose computer, special purpose computer, orother programmable data processing apparatus to produce a machine, suchthat the instructions, which execute via the processor of the computeror other programmable data processing apparatus, create means forimplementing the functions/acts specified in the flowchart and/or blockdiagram block or blocks. These computer readable program instructionsmay also be stored in a computer readable storage medium that can directa computer, a programmable data processing apparatus, and/or otherdevices to function in a particular manner, such that the computerreadable storage medium having instructions stored therein comprises anarticle of manufacture including instructions which implement aspects ofthe function/act specified in the flowchart and/or block diagram blockor blocks.

The computer readable program instructions may also be loaded onto acomputer, other programmable data processing apparatus, or other deviceto cause a series of operational steps to be performed on the computer,other programmable apparatus or other device to produce acomputer-implemented process, such that the instructions which executeon the computer, other programmable apparatus, or other device implementthe functions/acts specified in the flowchart and/or block diagram blockor blocks.

The flowchart and block diagrams in the Figures illustrate thearchitecture, functionality, and operation of possible implementationsof systems, methods, and computer program products according to variousembodiments of the present invention. In this regard, each block in theflowchart or block diagrams may represent a module, segment, or portionof instructions, which comprises one or more executable instructions forimplementing the specified logical function(s). It will also be notedthat each block of the block diagrams and/or flowchart illustration, andcombinations of blocks in the block diagrams and/or flowchartillustration, can be implemented by special purpose hardware-basedsystems that perform the specified functions or acts or carry outcombinations of special purpose hardware and computer instructions.

In the description and claims, each of the terms “substantially,”“essentially,” and forms thereof, when describing a numerical value,means up to a 20% deviation (namely, ±20%) from that value. Similarly,when such a term describes a numerical range, it means up to a 20%broader range-10% over that explicit range and 10% below it).

In the description, any given numerical range should be considered tohave specifically disclosed all the possible subranges as well asindividual numerical values within that range, such that each suchsubrange and individual numerical value constitutes an embodiment of theinvention. This applies regardless of the breadth of the range. Forexample, description of a range of integers from 1 to 6 should beconsidered to have specifically disclosed subranges such as from 1 to 3,from 1 to 4, from 1 to 5, from 2 to 4, from 2 to 6, from 3 to 6, etc.,as well as individual numbers within that range, for example, 1, 4, and6. Similarly, description of a range of fractions, for example from 0.6to 1.1, should be considered to have specifically disclosed subrangessuch as from 0.6 to 0.9, from 0.7 to 1.1, from 0.9 to 1, from 0.8 to0.9, from 0.6 to 1.1, from 1 to 1.1 etc., as well as individual numberswithin that range, for example 0.7, 1, and 1.1.

The descriptions of the various embodiments of the present inventionhave been presented for purposes of illustration, but are not intendedto be exhaustive or limited to the explicit descriptions. Manymodifications and variations will be apparent to those of ordinary skillin the art without departing from the scope and spirit of the describedembodiments. The terminology used herein was chosen to best explain theprinciples of the embodiments, the practical application or technicalimprovement over technologies found in the marketplace, or to enableothers of ordinary skill in the art to understand the embodimentsdisclosed herein.

In the description and claims of the application, each of the words“comprise,” “include,” and “have,” as well as forms thereof, are notnecessarily limited to members in a list with which the words may beassociated.

Where there are inconsistencies between the description and any documentincorporated by reference or otherwise relied upon, it is intended thatthe present description controls.

What is claimed is:
 1. A system comprising: at least one hardwareprocessor; and a non-transitory computer-readable storage medium havingstored thereon program instructions, the program instructions executableby the at least one hardware processor to: receive an input tensor,wherein said input tensor has a shape defined by [n₁, . . . , n_(k)],where k is equal to a number of dimensions that characterize the inputtensor, receive tile tensor metadata comprising at least: (i) a tiletensor shape defined by [t₁, . . . , t_(k)], and (ii) informationindicative of an interleaving stride to be applied with respect to eachdimension of said tile tensor, and construct an output tensor comprisinga plurality of said tile tensors, by applying a packing algorithm whichmaps each element of said input tensor to at least one slot location ofone of said plurality of tile tensors, based, at least in part, on saidtile tensor shape and said interleaving stride, wherein saidinterleaving stride results in non-contiguous mapping of said elementsof said input tensor, such that each of said tile tensors includes asubset of said elements of said input tensor which are spaced withinsaid input tensor according to said interleaving stride.
 2. The systemof claim 1, wherein said program instructions are further executable tostore said output tensor and said tile tensor metadata.
 3. The system ofclaim 1, wherein said program instructions are further executable tounpack said input tensor from said stored output tensor, based on saidtile tensor metadata.
 4. The system of claim 1, wherein said tile tensormetadata further comprises a replication parameter, and wherein saidpacking algorithm is configured to perform, based on said replicationparameters, replication of said elements of said input tensor, such thateach of said elements of said input tensor is mapped to multiple slotlocations along a dimension of one of said tile tensors.
 5. The systemof claim 1, wherein said program instructions are further executable to:receive a filter associated with a convolution computation over saidinput tensor; compute said convolution by applying a multiplyingoperator which multiplies said filter element-wise over each of saidtile tensors in said output tensor; apply a summation algorithm thatsums over the results of said multiplying; and output a result of saidapplying of said summation algorithm as a result of said convolution. 6.The system of claim 5, wherein said convolution is part of a neuralnetwork inference.
 7. The system of claim 1, wherein said tensor tilesare homomorphic encryption ciphertexts.
 8. A computer-implemented methodcomprising: receiving an input tensor, wherein said input tensor has ashape defined by [n₁, . . . , n_(k)], where k is equal to a number ofdimensions that characterize the input tensor; receiving tile tensormetadata comprising at least: (i) a tile tensor shape defined by [t₁, .. . , t_(k)], and (ii) information indicative of an interleaving strideto be applied with respect to each dimension of said tile tensor;constructing an output tensor comprising a plurality of said tiletensors, by applying a packing algorithm which maps each element of saidinput tensor to at least one slot location of one of said plurality oftile tensors, based, at least in part, on said tile tensor shape andsaid interleaving stride, wherein said interleaving stride results innon-contiguous mapping of said elements of said input tensor, such thateach of said tile tensors includes a subset of said elements of saidinput tensor which are spaced within said input tensor according to saidinterleaving stride.
 9. The computer-implemented method of claim 8,further comprising storing said output tensor and said tile tensormetadata.
 10. The computer-implemented method of claim 8, furthercomprising unpacking said input tensor from said stored output tensor,based on said tile tensor metadata.
 11. The computer-implemented methodof claim 8, wherein said tile tensor metadata further comprises areplication parameter, and wherein said packing algorithm is configuredto perform, based on said replication parameters, replication of saidelements of said input tensor, such that each of said elements of saidinput tensor is mapped to multiple slot locations along a dimension ofone of said tile tensors.
 12. The computer-implemented method of claim8, further comprising: receiving a filter associated with a convolutioncomputation over said input tensor; computing said convolution byapplying a multiplying operator which multiplies said filterelement-wise over each of said tile tensors in said output tensor;applying a summation algorithm that sums over the results of saidmultiplying; and outputting a result of said applying of said summationalgorithm as a result of said convolution.
 13. The computer-implementedmethod of claim 12, wherein said convolution is part of a neural networkinference.
 14. The computer-implemented method of claim 8, wherein saidtensor tiles are homomorphic encryption ciphertexts.
 15. A computerprogram product comprising a non-transitory computer-readable storagemedium having program instructions embodied therewith, the programinstructions executable by at least one hardware processor to: receivean input tensor, wherein said input tensor has a shape defined by [n₁, .. . , n_(k)], where k is equal to a number of dimensions thatcharacterize the input tensor; receive tile tensor metadata comprisingat least: (i) a tile tensor shape defined by [t₁, . . . , t_(k)], and(ii) information indicative of an interleaving stride to be applied withrespect to each dimension of said tile tensor; and construct an outputtensor comprising a plurality of said tile tensors, by applying apacking algorithm which maps each element of said input tensor to atleast one slot location of one of said plurality of tile tensors, based,at least in part, on said tile tensor shape and said interleavingstride, wherein said interleaving stride results in non-contiguousmapping of said elements of said input tensor, such that each of saidtile tensors includes a subset of said elements of said input tensorwhich are spaced within said input tensor according to said interleavingstride.
 16. The computer program product of claim 15, wherein saidprogram instructions are further executable to store said output tensorand said tile tensor metadata.
 17. The computer program product of claim15, wherein said program instructions are further executable to unpacksaid input tensor from said stored output tensor, based on said tiletensor metadata.
 18. The computer program product of claim 15, whereinsaid tile tensor metadata further comprises a replication parameter, andwherein said packing algorithm is configured to perform, based on saidreplication parameters, replication of said elements of said inputtensor, such that each of said elements of said input tensor is mappedto multiple slot locations along a dimension of one of said tiletensors.
 19. The computer program product of claim 15, wherein saidprogram instructions are further executable to: receive a filterassociated with a convolution computation over said input tensor;compute said convolution by applying a multiplying operator whichmultiplies said filter element-wise over each of said tile tensors insaid output tensor; apply a summation algorithm that sums over theresults of said multiplying; and output a result of said applying ofsaid summation algorithm as a result of said convolution.
 20. Thecomputer program product of claim 19, wherein said convolution is partof a neural network inference.